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ABSTRACT 


A  weather  radar  signal  simulator  that  produces  an  output  consisting  of  a  vector  of 
I  and  Q  values  representing  the  radar  return  permits  investigation  of  the  perfonnance  of 
different  estimators  for  the  weather  signal  parameters  and  their  sensitivity  when  varying 
radar  parameters  and  precipitation  models.  Although  several  empirical  statistical  models 
are  available  to  describe  precipitation  behavior,  the  creation  of  a  physical  model  enables 
adaptation  to  actual  data  (e.g.  rain  rate,  wind  shears)  thereby  making  it  possible  to  apply 
and  examine  different  scanning  schemes,  especially  rapid  scanning  schemes.  A  physical 
model  allows  gradual  improvements  to  realism  to  study  the  effects  on  the  radar  return  for 
different  phenomena.  A  Weather  Radar  Signal  Simulator  has  been  developed  in 
MATLAB.  Several  different  functionalities  have  been  implemented  allowing  for  stepped 
frequency,  multiple  PRFs,  pulse  compression  using  a  chirp,  and  variation  of  both  weather 
and  radar  input  parameters.  Post  processing  capabilities  include  autocorrelation  and  FFT 
(for  single  PRF  only);  estimation  of  weather  parameters  such  as  reflectivity  factor,  Z; 
average  doppler,  radial  velocity,  and  velocity  spread;  pedagogical  plots  including  a 
Phasor  plot  of  phase  change  over  time  and  a  velocity  histogram,  instantaneous  observed 
reflectivity  and  power  for  each  pulse  over  time. 
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EXECUTIVE  SUMMARY 


Electromagnetic  emission  is  increasing  in  all  surroundings  not  only  in  the  military 
arena  making  it  necessary  to  understand,  interpret,  and  sometimes  filter  information 
created  by  man  or  nature.  In  an  attempt  to  learn  more  in  one  research  area,  studying 
other  areas  might  provide  the  insight  needed  to  break  new  ground.  Electromagnetic 
returns  from  weather  are  in  one  area  of  radar  science  regarded  as  clutter  while  in  another 
it  as  a  way  of  gaining  information  and  making  estimates  and  predictions.  Exploring  the 
weather  side  of  radar  science  provides  an  insight  into  weather  returns  as  source  of 
information  and  also  improves  knowledge  about  effects  that  can  lead  to  a  better 
understanding  of  how  to  overcome  unwanted  weather  returns. 

Creating  a  weather  radar  signal  simulator  that  produces  an  output  consisting  of  a 
vector  of  I  and  Q  values  representing  the  radar  return,  pennits  investigation  of  the 
performance  of  different  estimators  for  the  weather  signal  parameters  and  their  sensitivity 
when  varying  radar  parameters  and  precipitation  models.  In  weather  radar  applications 
the  return  signal  samples  can  be  used  to  estimate  reflectivity,  velocity  and  velocity 
spread.  These  are  the  zeroth,  first  and  second  moments  of  the  Doppler  spectrum. 
Reflectivity  can  be  used  to  estimate  the  rain  rate  of  the  specific  volume  return.  A  number 
of  mathematical  models  relating  rain  rate  to  reflectivity  exist  and  can  be  used  when  an 
estimated  power  return  is  available.  There  are  also  a  number  of  measures  and  statistics 
for  precipitation  that  have  resulted  in  functions  describing  rain  drop  size  distribution  for  a 
certain  rain  rate,  terminal  velocities  of  rain  drops,  etc. 

Limiting  the  amount  of  samples  required  to  accomplish  small  variance  when 
estimating  average  power  calls  for  independent  samples  meaning  that  the  return  samples 
are  de-correlated.  On  the  other  hand,  measuring  Doppler  calls  for  coherent,  correlated 
samples,  this  typically  requires  a  higher  pulse  repetition  frequency  (PRF).  The  velocity 
spread  effects  the  weather  signal  correlation  time  and  therefore  drives  sample  time 
required  to  arrive  at  an  estimate  with  specific  variance.  Although  several  empirical 
statistical  models  are  available  to  describe  precipitation  behavior,  the  creation  of  a 
physical  model  enables  adaptation  to  actual  data  (e.g.  rain  rate,  wind  shears)  thereby 

xvii 


making  it  possible  to  apply  and  examine  different  scanning  schemes,  especially  rapid 
scanning  schemes,  as  well  as  to  create  output  graphical  representations  that  can  serve  a 
pedagogical  purpose. 

This  research  was  undertaken  to  build  a  weather  radar  signal  simulator  delivering 
I  and  Q  values  from  a  physical  representation  of  rain  drops,  which  has  been 
accomplished.  A  previously  developed  simulator  was  initially  evaluated  and  verified 
serving  as  the  basic  building  block  when  developing  the  final  version.  The  Radiation 
Integrals  were  used  to  derive  the  electromagnetic  scattering  from  rain  drops 
approximated  as  dielectric  spheres.  Without  accounting  for  multiple  scattering, 
attenuation,  or  coupling  effects,  the  scattered  electric  field  from  each  drop  is  summed  for 
every  pulse  providing  the  RCS  of  the  rain  in  the  radar  resolution  cell  as  well  as  a  complex 
voltage  return  representation.  Using  the  Marshall-Palmer  drop  size  distribution,  an  initial 
position  for  every  drop  is  determined  within  an  extended  radar  resolution  volume 
enabling  the  drops  to  fall  through  the  radar  beam  for  all  samples.  Between  pulses  the 
drops  are  moved  according  to  their  respective  terminal  velocities  based  on  the  Atlas- 
Ulbrich  approximation  and  wind.  To  allow  for  width  in  the  velocity  spectrum,  a 
Gaussian  spread  is  applied  to  drop  velocities  that  determine  how  the  average  velocity  is 
spread.  Also  a  Gaussian  spread,  that  is  size  dependent,  is  applied  to  the  velocity  and  it 
represents  the  spread  of  velocity  due  to  difference  in  size.  The  simulator  power  return  is 
weighted  using  the  Radar  Range  Equation  for  point  targets,  considering  all  drops  as  point 
targets  within  the  resolution  cell.  The  zeroth  moment  estimate  is  derived  by  averaging  N 
power  samples. 

Several  different  functionalities  have  been  implemented  allowing  for  stepped 
frequencies,  multiple  PRFs,  pulse  compression  using  a  chirp,  and  varying  input 
parameters.  Post  processing  capabilities  include  autocorrelation  and  FFT  (only  for  single 
PRF);  evaluation  of  weather  parameter  estimators  such  as  average  reflectivity  factor,  Z; 
average  doppler,  radial  velocity,  and  velocity  spread.  Pedagogical  plots  including  a 
Phasor  plot  of  phase  change  over  time  and  a  velocity  histogram,  instantaneous  observed 
reflectivity  and  power  for  each  pulse  over  time,  have  also  been  implemented. 


I.  INTRODUCTION 


A.  BACKGROUND 

Electromagnetic  emission  is  increasing  in  all  surroundings  not  only  in  the  military 
arena  making  it  necessary  to  understand,  interpret,  and  sometimes  filter  information 
created  by  man  or  nature.  In  an  attempt  to  learn  more  in  one  research  area,  studying 
other  areas  might  provide  the  insight  needed  to  break  new  ground.  Electromagnetic 
returns  from  weather  are  in  one  area  of  radar  science  regarded  as  clutter  while  others  use 
it  as  a  way  of  gaining  information  and  make  estimates  and  predictions.  Exploring  the 
weather  side  of  radar  science  will  not  only  provide  an  insight  into  weather  returns  as 
source  of  information  but  also  improve  knowledge  about  effects  that  can  lead  to  a  better 
understanding  of  how  to  overcome  unwanted  weather  returns. 

Creating  a  weather  radar  signal  simulator  that  produces  an  output  consisting  of  a 
vector  of  I  and  Q  values  representing  the  radar  return,  will  permit  investigation  of  the 
performance  of  different  estimators  for  the  weather  signal  parameters  and  their  sensitivity 
when  varying  radar  parameters  and  precipitation  models.  In  weather  radar  applications 
the  return  signal  samples  can  be  used  to  estimate  reflectivity,  velocity  and  velocity 
spread.  These  are  the  zeroth,  first  and  second  moments  of  the  Doppler  spectrum. 
Reflectivity  can  be  used  to  estimate  the  rain  rate  of  the  specific  volume  return.  A  number 
of  mathematical  models  relating  rain  rate  to  reflectivity  exist  and  can  be  used  when  an 
estimated  power  return  is  available.  There  are  also  a  number  of  measures  and  statistics 
for  precipitation  that  have  resulted  in  functions  describing  rain  drop  size  distribution  for  a 
certain  rain  rate,  terminal  velocities  of  rain  drops,  etc. 

To  minimize  the  number  of  samples  required  to  realize  small  variance  when 
estimating  average  power  independent  samples  are  required  meaning  that  the  return 
samples  are  de-correlated.  On  the  other  hand,  measuring  Doppler  calls  for  coherent, 
correlated  samples,  this  typically  requires  a  higher  pulse  repetition  frequency  (PRF).  The 
velocity  spread  effects  the  weather  signal  correlation  time  and  therefore  drives  sample 
rate.  Although  several  empirical  statistical  models  are  available  to  describe  precipitation 
behavior,  the  creation  of  a  physical  model  enables  adaptation  to  actual  data  (e.g.  rain  rate, 
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wind  shears)  thereby  making  it  possible  to  apply  and  examine  different  scanning 
schemes,  especially  rapid  scanning  schemes,  as  well  as  to  create  output  graphical 
representations  that  can  serve  a  pedagogical  purpose.  The  output  signal  should  therefore 
correspond  to  the  values  expected  using  the  empirical  models.  The  simulator  should  also 
enable  input  of  data  to  create  variation  in  reflectivity,  velocity  spread  and  Doppler. 

The  simulator  is  based  on  the  reflection  of  each  drop  creating  a  scattered  electric 
field.  A  Marshall-Palmer  drop  size  distribution  is  used.  The  sum  of  the  fields  scattered 
by  all  drops  at  any  particular  instant  of  time  gives  the  instantaneous  scattered  field,  total 
radar  cross-section  and  instantaneous  received  power.  Motion  is  imparted  to  each  drop 
using  a  tenninal  velocity  based  on  drop  size  and  wind  field.  The  motion  leads  to  a 
rearrangement  of  drops  in  the  sample  volume  and  thus  to  a  new  value  of  instantaneous 
received  weather  signal  power  at  the  next  sample  time.  The  effect  of  electromagnetic 
coupling  among  drops  is  investigated. 


B.  OBJECTIVE 

The  main  objective  was  to  build  a  weather  radar  signal  simulator  that  produces  an 
output  consisting  of  a  vector  of  I  and  Q  values  of  the  radar  return.  From  this,  average 
power,  Doppler  and  Doppler  spread  can  be  derived.  To  guide  the  research  and  provide 
tools  for  verification  and  validation  of  the  simulator  the  following  questions  were 
addressed: 

1 .  Is  it  possible  to  build  a  weather  radar  signal  simulator  based  on  a  physical 
model  of  a  spatial  region  containing  raindrops? 

2.  Can  realistic  motion  that  depends  upon  the  type  of  weather  system  be 
imparted  to  the  raindrops? 

3.  Can  the  output  of  the  weather  radar  simulator  be  used  to  study  spatial  and 
temporal  sampling  schemes  for  constant  or  stepped  frequency  sampling 
pulses? 
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4.  Can  the  weather  radar  simulator  be  used  to  study  the  use  of  pulse 
compression  and  range  averaging  as  a  means  of  rapidly  obtaining 
independent  samples? 

5.  Can  the  output  of  the  weather  radar  signal  simulator  be  used  to  produce 
displays  of  pedagogical  interest? 

6.  Can  the  output  of  the  weather  radar  signal  simulator  be  used  to  study  the 
performance  of  estimators  for  the  first  three  Doppler  moments? 

7.  Can  the  output  of  the  weather  radar  signal  simulator  be  used  to  study  the 
utility  of  estimators  for  higher  order  Doppler  moments? 

C.  APPROACH 

The  simulator  was  developed  by  building  modules  that  represented  a  functionality 
or  physical  behavior,  testing  and  verifying  the  module,  and  then  placing  the  module  into 
the  final  version  of  the  simulator.  Having  a  prototype  from  France  provided  valuable 
ideas  and  experience  for  developing  the  simulator.  Initially  this  version  was  translated, 
evaluated,  and  respective  function  or  module  output  was  verified  against  empirical 
models  based  on  statistics  of  real  life  measurements.  When  all  improvements  were 
implemented  to  the  French  version,  new  modules  were  created  aiming  at  addressing  the 
research  questions  posed.  The  final  product  was  optimized  to  give  accurate  results 
structured  and  presented  in  a  pedagogical  manner.  During  the  development  Graphical 
User  Interfaces  (GUIs)  were  of  secondary  concern. 

When  the  simulator  was  working,  i.e.  it  has  passed  the  necessary  testing  for 
verification;  the  application  phase  began.  Research  questions  three  to  seven  were  then 
addressed.  Results  were  again  compared  with  empirical  models  built  on  statistics  of  real 
life  measurements.  Questions  not  addressed  were  fonnulated  for  future  investigation  by 
others. 

Although  the  subject  title  implies  that  meteorology  will  be  extensively  covered, 
the  goal  of  this  work  was  to  provide  a  tool  to  evaluate  the  effects  on  the  radar  output  of 
different  weather  systems.  Precipitation  physics  was  used  to  build  a  physical  model,  but 
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initially  weather  system  input  was  included  only  as  a  conceptual  mathematical  vector 
model  to  provide  flexibility  to  later  implement  real  models  of  weather. 


D.  RELATED  WORK 

M.  Gosset,  J.  Nicol,  and  A.  Sanchez  submitted  an  abstract  for  The  6th  Conference 
on  Hydrological  Applications  of  Weather  Radar  describing  a  simulator  with  similar  focus 
[1].  After  personal  contact  with  M.  Gosset  we  learned  that  the  simulator  was  unfinished 
and  later  the  MATLAB  code  was  sent  to  us  to  develop  at  NPS  in  cooperation  with  M. 
Gosset. 

J.S.  Marshall  and  Walter  Hitschfeld  [2]  thoroughly  investigated  theoretical  effects 
such  as  fluctuation  echoes  from  randomly  distributed  scatterers  providing  mathematical 
proofs  and  suggestions  for  pulsing  and  scanning  schemes  to  improve  estimate  accuracy. 
The  paper  defines  probability  distributions  useful  to  interpret  fluctuating  returns  and 
describes  different  methods  to  achieve  independent  data  using  difference  in  time,  range, 
and  frequency,  to  mention  a  few. 

Richard  J.  Doviak  and  Dusan  S.  Zmic  [3]  provide  extensive  guidance  in  the 
weather  radar  area  presenting  tools  and  explanations  of  all  aspects  investigated  in  this 
thesis.  Most  background  data  has  originated  from  thoughts,  theories,  and  descriptions 
provided  by  their  fantastic  book.  Building  a  simulator  calls  for  in  depth  interpretation  of 
radar  hardware  solutions,  which  in  large  part  can  be  found  in  their  work. 

Finally,  John  B  Sandifer  [4]  developed  scanning  strategies  for  research  and 
operational  applications.  Interesting  temporal  schemes  were  presented  aiming  to  provide 
methods  to  rapidly  scan  and  update  large  volumes  using  Phased  Array  Doppler  radars. 
The  use  of  stepped  frequency  to  acquire  independent  samples  without  having  to  wait  for 
the  weather  to  reshuffle  enough  for  de-correlation  is  explored  and  presented  as  one 
successful  way  of  achieving  a  fast  scanning  capability. 
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E.  THESIS  OUTLINE 

The  thesis  is  organized  as  follows: 

Chapter  II,  titled  Weather  Radar  and  Weather/Rain  theory,  introduces  the  basic 
principles  of  radar  and  the  specific  functions  in  weather  applications.  The  physics  of 
precipitation  is  summarized  with  emphasis  on  electromagnetic  effects  (scatter),  motion 
behavior,  and  drop  size  distribution.  The  chapter  also  describes  Radar  Cross  Section 
(RCS)  of  rain.  The  basis  of  the  simulator  is  the  radar  signal  return  from  the  rain,  which  is 
why  the  electromagnetic  effects  will  be  examined  to  provide  the  foundation  for  the 
simulator  design. 

Chapter  III  covers  simulator  design  considerations  and  the  building  of  the 
simulator.  In  this  chapter  the  scope  and  the  limitations  of  the  simulator  are  presented 
based  upon  the  research  questions  that  need  to  be  addressed  and  the  abilities  of  the 
simulator  building  tool  of  choice  (MATLAB).  The  chapter  further  describes  the  building 
of  the  simulator.  This  chapter  covers  both  the  development  and  testing  of  the  simulator. 
MATLAB  is  the  computer  language  of  choice  and  the  simulator  is  initially  developed  in 
modules,  to  simplify  development  and  testing.  Ending  the  chapter  the  simulator  parts  are 
described. 

Chapter  IV  provides  analysis  of  applications.  The  final  test  or  validation  of  the 
simulator  where  the  research  questions  3-7  are  addressed. 

Chapter  V  presents  the  results  from  both  the  development  and  the  applications 
covered  in  the  thesis. 


Chapter  VI  contains  conclusions  and  recommendations  for  future  work 

Appendix  A  Derivation  of  RCS  for  single  drop 

Appendix  B  Microwave  Studio  simulator  setup 

Appendix  C  Simulator  setup  for  the  Weather  Radar  Simulator 

Appendix  D  MATLAB  code 
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II.  WEATHER  RADAR  AND  WEATHER/RAIN  THEORY 


Laying  the  foundation  for  the  simulator,  the  basic  theory  is  covered  in  this  chapter 
ending  with  a  statement  of  ground  truth  for  later  comparison  and  verification. 


A.  RADAR  PRINCIPLES 

Radar  or  RAdio  Detecting  And  Ranging  uses  the  return  of  transmitted 
electromagnetic  energy  from  a  transmitter  to  detect  and  locate  objects  of  interest. 
Although  the  basic  principles  of  radar  have  been  known  for  decades  new  applications  are 
developed  constantly  as  a  consequence  of  the  exploration  of  different  areas  of  interest. 
The  basic  function  of  a  radar  is  provided  by  a  transmitter  that  generates  an 
electromagnetic  signal  into  space  by  the  means  of  an  antenna.  Objects  in  the  signal  path 
scatter  or  reflect  part  of  the  power  projected  by  the  transmitter  back  towards  the  receiver. 
By  measuring  the  return  in  terms  of  angle  of  incidence  and  time  to  the  return,  power 
levels,  etc.  different  factors  can  be  estimated  depending  on  the  application.  [5] 

To  be  able  to  detect  movement  or  to  calculate  average  power,  more  than  one 
return  pulse  must  be  measured  leading  to  the  idea  of  Pulse  Repetition  Frequency  (PRF). 
By  selecting  different  PRFs  some  quantities  will  be  become  ambiguous  while  others 
become  the  unambiguous. 

1.  Estimation  of  Range 

Finding  the  range  to  the  object  of  interest  is  merely  a  question  of  measuring  the 
time  it  takes  for  the  signal  to  travel  from  the  transmitter  to  the  target  and  back  to  the 
receiver,  or  more  precisely  from  and  to  the  antenna.  Knowing  that  the  electromagnetic 
energy  travels  at  the  speed  of  light,  c,  which  in  free  space  is  approximately  3x10s  m/s, 
the  distance  R  is  given  by 

Rmf  p-‘> 

where  T  is  the  time  from  the  antenna  to  the  object  and  back  to  the  antenna. 
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To  ensure  unambiguous  measurements  of  range  the  time  difference  between  two 
consecutive  pulses  must  be  at  least  the  time  it  takes  for  the  signal  to  travel  to  the  range  of 
interest  and  back.  The  maximum  ambiguous  range,  Ru  is  thereby 


r  = 


2  L 


where  Trep  is  the  pulse  repetition  period  and  fp  the  PRF. 


(2.2) 


2  Estimation  of  Power  Density 

a.  Antenna  Gain  and  Beam  Shapes 

An  important  part  of  a  radar  system  is  the  antenna.  The  antenna  is  the 
impedance  matcher  making  it  possible  to  transmit  and  direct  energy  into  space  with 
minimal  losses.  Some  antennas  provide  the  ability  to  direct  the  energy,  making  the 
power  density  higher  in  a  limited  solid  angle  while  being  lower  elsewhere.  This  ability, 
called  directivity,  is  converted  to  a  measure  of  power  per  unit  solid  angle  radiated  in  a 
particular  direction  and  is  compared  to  the  same  measure  but  without  directivity.  The 
ratio  is  called  the  directive  gain  of  an  antenna  and  is  a  function  of  direction  and  it  does 
not  include  losses  in  the  antenna. 

Apart  from  resulting  in  a  higher  Gain,  the  ability  to  shape  the  beam  can  be 
helpful  to  enhance  angular  resolution,  since  it  will  limit  the  volume  of  illumination.  An 
isotropic  antenna  has  a  power  density  at  a  range  R  of 


p=-A 


\nR~ 


(2.3) 


where  Pt  is  the  transmitted  power,  while  an  antenna  with  gain  has  the  power  density  of 


4  7rR 


2  * 


(2.4) 


b.  Radar  Range  Equation 

To  be  able  to  estimate  the  return  power  from  an  object,  a  parameter  called 
Radar  Cross  Section  (RCS)  must  be  defined.  The  RCS  can  be  defined  as  [6] 
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_  Power  reflected  to  receiver  per  unit  solid  angle 
Incident  power  density/47i 

The  scattered  density  power  at  the  receiving  antenna  is  given  by 

P,G  < 


P  =— •- 

ant 


[W/m  ]. 


4  X  R1  4k  R2 

Accounting  for  the  capture  area  of  the  radar  antenna,  Aeff  ,  results  in  received  power 

By  using  the  relationship  between  transmit  gain  and  receive  effective  area, 

4  kA„ 


(2.5) 


(2.6) 


G  = 


A2 


(2.7) 


where  A  is  the  radar  wavelength,  received  power  can  be  expressed  as 

_  PtGtGraXr 

rec  /  \3  n4  (2-0) 

(4 k)  R 4 

ignoring  all  system  and  propagation  losses.  Since  losses  are  expected  let  us  for  now  add  a 
one-way  propagation  loss  term  L  rendering 


p  _  PfrG'XTA2 


(4  kYr4L2 


(2.9) 


3.  Estimation  of  Velocity 

Pulsed  radar  can,  by  measuring  the  phase  change  between  consecutive  pulses, 
estimate  radial  velocity.  The  effect,  called  Doppler,  is  the  frequency  change  that  occurs 
due  to  the  relative  movement  of  a  target  with  respect  to  the  radar.  Considering  a  two-way 
propagation  path  renders  a  total  phase  change  of 

2  R 

(/)  =  2k - .  (2.10) 

A 

In  terms  of  angular  frequency  the  result  is 
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(2.11) 


d(b  An  dR  Anvr  .  . 

CU  =  —  = - = - L  =  In  f, 

dt  A  dt  A 

where  vr  is  the  radial  velocity  and  fd  is  the  doppler  frequency  [4]. 


From  a  pulse  pair,  using  a  coherent  detector,  the  phase  change  can  be  extracted 


using 


V  =1  +  jO  =  A  ej 

n  n  J  n 

arg  [K,}  =  0n 


v  V*  =  A  A  eJ 

v  n+V  n  ^n+X^n^ 


arg{F„+iF;}  =  ^+i-^  =  ^ 

leading  to  the  single  pulse  pair  velocity  estimate 


(2.12) 


A  Sty  A  (-  *1  Af  (■  *1 

^  = - -  = - arg{F„+1F„  }  =  ^arg [Vn+lVn  } 


AnSt  AnT„ 


(2.13) 


where  T  is  time  between  pulses  and  /  is  pulse  repetition  frequency. 


B.  WEATHER  RADAR 

The  radar  application  that  will  be  explored  in  this  research  is  the  ability  to  use  the 
radar  to  estimate  weather  parameters  such  as  rain  reflectivity  and  wind  speed.  Most  radar 
applications  consider  rain  as  clutter  and  find  ways  to  filter  out  those  effects.  In  weather 
applications  the  scatter  from  precipitation  is  used  to  compute  the  estimates  mentioned 
earlier. 


1.  Sample  Correlation 

The  return  signal  power  level  from  weather  will  vary  over  time.  To 
compute  an  estimate  of  reflectivity  with  small  variance  a  number  of  independent  samples 
are  summed  and  averaged  which  means  that  the  target  volume  must  reshuffle  enough  to 
ensure  de-correlation  during  the  observation  time.  The  power  spectral  density  of  a 
meteorological  signal  is  approximately  Gaussian  [7]  and  can  be  written  as 
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-fczfi 

2  cl 

S(/)  =  V'  1  (2.14) 

where  /  is  the  doppler  frequency  and  <72  is  the  doppler  spectrum  variance.  Taking  the 
Fourier  transform  of  (2.14)  and  nonnalizing  results  in  the  correlation  coefficient 


p(T) 


T~ 

2(7 ] 


(2.15) 


where  o]  is  the  time  spectrum  variance  and  is  related  to  o2  by 


and 


1 

2k<j  f 


(2.16) 


(2.17) 


By  using  the  correlation  coefficient  a  measure  of  minimum  difference  in  time  can 
be  estimated  ensuring  dependent  or  independent  samples.  Nathanson  claims  that  for 
independent  sampling  /l(r)  <  0.02  is  required  whereas  /l(r)  >  0. 1 5  will  insure 
dependence  [7].  Doviak  and  Zmic  [3]  set  a  correlation  threshold  for  coherence  at 
p1  (t)  >  e  1  which  corresponds  to  a  value  of  p{j)  >  0.6  .  Using  this  threshold  renders 


Ts  <0.25—  (2.18) 

nov 


for  highly  correlated  samples.  For  independent  samples  using  [7]  T  >0.70 -  is 

required. 


2.  Depth  of  Volume 

Estimating  range  to  a  volume  target,  compared  to  a  point  target,  still  involves 
measuring  travel  time.  The  range  to  every  raindrop  will  not  be  detennined  since  the 
return  signal  will  be  the  sum  of  several  scatterers.  Measuring  a  volume  takes  pulsewidth 
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into  account  as  one  of  the  parameters  detennining  its  size.  A  pulsewidth  of  r  renders  a 
depth  of  the  return  volume  of 


(ly 


CT 

T' 


(2.19) 


3.  Estimation  of  Power  Density  and  Reflectivity  Z 

a.  Antenna  Gain  and  Beam  Effective  Solid  Angle 

The  directivity  of  an  antenna  is  a  function  of  both  azimuth  and  elevation 
angles.  The  radar  beam  has  a  width  that  has  to  be  specified  to  enable  an  estimation  of  the 
magnitude  of  the  volume  return.  Estimates  of  a  point  target  involve  the  angle  off  bore 
sight  to  get  the  accurate  gain  level,  while  a  volume  target,  like  rain,  takes  into  account  the 
total  gain  over  the  volume  of  interest.  Estimation  of  the  gain  for  a  volume  target  uses  the 
antenna  pattern  function  f(6,(p )  integrated  over  the  solid  angle.  Since  the  radar 

equation  is  a  power  relation  the  pattern  function  needs  to  be  squared,  also,  due  to  the  two 
way  propagation,  the  power  pattern  function  is  squared  again  leading  to 

j|/(0,^)|4ja  (2.20) 

4  n 


Using  a  Gaussian  Pattern  Model  provides  a  good  approximation  of  the 
main  beam  between  the  3  dB  points.  Extending  it  to  both  azimuth  and  elevation  renders 

_( e2 1  j 

\f(e,4>f=e^  A  |0|<U„  (2.21) 

Q‘ -  rff 

where  y1  =  — 1 —  and  S1  =  — — ,  and  0X ,  <p}  are  the  half  power  beamwiths  in  the 
4  In  2  4  In  2 

orthogonal  planes.  Performing  the  integration  over  the  solid  angle  renders  [8] 


'  w1 

2(U 

v  r2 

s2  j 

d6d(j)  ■ 


4  K 


nOx(px 
8  In  2 


(W) 


(2.22) 


where  E,  <  0.034  and  represents  the  contribution  from  the  sidelobes.  It  corresponds  to 
0.15  dB  error  allowing  for  an  approximation  ignoring  the  term  {\-f)  rendering 


12 


J|/(M 


4  _KOf>x 

_  8  In  2 


(2.23) 


b.  Weather  Radar  Range  Equation  and  Estimation  of  Reflectivity 
To  define  the  weather  radar  range  equation  we  return  to  the  original 


formula 


_  PtGtGrai l2 
rec  (4 n)3R4  ' 


(2.24) 


In  the  weather  radar  equation  RCS  (a)  is  replaced  by  RCS  density  // ,  which  is  the 
expected  radar  cross  section  density  per  unit  volume.  Since  the  target  is  a  volume  target, 
the  expected  return  needs  to  be  integrated  over  the  whole  volume,  which  can  be  divided 
into  the  solid  angle  and  depth  of  the  return.  Letting  the  power  level  at  the  target  be 
PG  PG 

R  „„  = — ’—fdA  =  -J—LdQ  and  integrating  the  return  power  over  the  return  volume 
target  AnR2  An 

renders 


R+— 

2 


P  _  f  jP  _  ^2  G,Gr  r-  dR  n  ,  ,u4  _  PtA~GtGrcr0^ 

p~ -\dp™— (tof-  J  da-  I024»j/ein2  n  (2'25) 


recognizing  that  R  »  — . 


Radar  meteorologists  use  reflectivity  factor,  Z,  instead  of  instead  of  RCS 


density  where 


n=Y  Kfz' 


(2.26) 


Z  e— Vp6. 

Substituting  (2.26)  in  (2.25)  and  solving  for  Z  renders  the  reflectivity  estimator 


(2.27) 


1 024^7? 2  In  2  Z4  _  -  1024R2ln2Z2 

PtXrGtGrcTdf)x  n5  Kw  2  PfjfjctO^n1  Kw  2 


(2.28) 
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For  an  assumed  continuous  distribution  of  drops  sizes  with  diameter  density 
MID)  (2.27)  becomes 

Z  =  |"m“  N(D)D6dD.  (2.29) 

Solving  the  integral  form  for  an  assumed  Marshall-Palmer  Drop  Size  Distribution, 
DSD  (Equation  (2.35)),  yields 

Z  =  }l?N(D)dD  =  6W°  (2.30) 

o  (4.  ITT0'21) 


Drops  do  not  exist  in  sizes  ranging  up  to  infinite  diameter  (see  later  sections  for 
details)  but  are  likely  to  be  limited  by  a  maximum  diameter.  Doviak  and  Zrnic  [3]  state 
that  Equation  (2.30)  overestimates  Z  and  suggest  a  truncated  DSD  with  D  =  Dmax .  In  this 
case 


Z  = 


Nn 


(4.  ITT0'21) 


T(7,o) 


(2.31) 


where  a  =  ADmax  and  y{l,a)  is  the  incomplete  Gamma  function.  Equations  (2.30)  and 

(2.31)  provide  an  analytical  means  for  computing  Z  and  providing  a  reference  as  ground 
truth  reflectivity  against  which  simulator  estimates  can  be  compared. 


4.  Estimation  of  Velocity 

Measuring  phase-change  over  time  for  a  volume  target  demands  coherence 
between  pulse  pairs.  Phase  change  will  also  occur  due  to  other  effects  rather  than  radial 
movement.  To  obtain  an  estimate  of  radial  velocity  with  small  variance,  a  number  of 
sample  pairs  need  to  be  summed  and  averaged.  Building  on  (2.13)  we  obtain  the  multiple 
pulse  pair  estimator 


1  N-\ 

— Yi 

n- itr 


n+l 


N-\  An 


Zar§R+i  v*,\ 


(2.32) 
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5.  Estimation  of  Velocity  Spread 

An  estimate  of  the  velocity  spread  can  be  computed  using  the  velocity  samples 
obtained  from  a  pulse  pair  comparison  (Equation  (2.13)).  The  estimator  for  sample 
variance  is  given  by 


and  the  estimator  for  standard  deviation  is 


(2.33) 


(2.34) 


C.  PHYSICS  OF  PRECIPITATION 

Why  and  how  raindrops  fonn  and  what  constitutes  their  shape  and  behavior  is 
beyond  the  scope  of  this  thesis,  however,  some  basics  in  precipitation  physics  are  needed 
to  develop  a  physical  representation  of  rain  for  the  simulator. 

1.  Size  and  Shape 

All  drops  examined  are  assumed  to  be  falling  at  or  close  to  their  terminal  velocity. 
Small  drops,  with  a  diameter  D  <  0.35  mm,  are  spherical  in  shape.  Drops  ranging  from 
0.35  <  D  <4  mm  have  progressively  flattened  bases,  and  can  be  approximated  by 
spheroids  [2].  Larger  raindrops  tend  to  break  up  after  collision  and  drops  larger  than 
about  10  mm  in  diameter  are  unstable  and  break  up  even  without  collision  [9].  Also 
drops  can  cluster  which  will  affect  the  drop  size  distribution.  This  phenomenon  will  not 
be  covered  in  this  thesis.  For  further  reading  [10]  is  suggested. 

2.  Drop  Size  Distribution,  DSD 

RCS  density  is  the  average  radar  cross  section  density  per  unit  volume,  as 
mentioned  previously,  and  for  independent  scattering  is  the  sum  of  the  RCS  of  all 
scatterers  in  the  return  volume  divided  by  the  volume.  The  drop  size  distribution 
represents  the  diameter  density  (mm/m  )  of  drops  for  all  diameters.  The  Marshall-Palmer 
drop  size  distribution  [11],  derived  empirically,  provides  a  general  model, 
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(2.35) 


TV (D)  =  N0  exp(-AD), 
A  =  4AR~°2lmm~\ 

N0  =8x1  O'  in  mm  1 . 


Applying  truncation  as  suggested  in  [3]  with  a  maximum  drop  diameter  of  Dmax 

yields 


N0  exp(-AD), 
0, 


D  <D 

max 

D . <  D 


(2.36) 


3.  Terminal  Velocities 

The  terminal  velocity  of  a  water  drop  is  dependent  on  the  size  of  the  drop.  Atlas 
and  Ulbrich  [12]  derived  an  expression  that  relates  terminal  velocity  to  the  diameter  of 
the  water  drop  as 

v(Z>)  =  1767 D°J7\cm/s\  «  386.6 D°f[m/s]  (2.37) 

which  is  valid  for  5xlCT4  <  Z)  <  5xlCT3,  [m] .  The  result  has  been  verified  against 
measured  data  [3],  Fall  speed  is  also  dependent  on  air  pressure  but  that  is  a  second  order 
effect  not  considered  here. 

4.  Velocity  Spectrum  Width 

There  are  several  mechanisms  contributing  to  the  spread  of  the  velocity  spectrum 
both  from  nature  and  the  radar  itself.  The  radar  antenna  motion,  turbulence,  differences 
in  fall  speed  due  to  drop  size  and  wind  shears  all  contribute  to  the  spread  of  the  radial 
velocity  of  precipitation.  Physically  modeling  all  contributors  demands  great  insight  into 
both  behavior  and  effects,  and  is  beyond  the  scope  of  this  thesis.  However,  to  control  de- 
correlation  time  a  spread  must  be  introduced  in  the  simulator.  This  is  accomplished  by 
imparting  a  random  Gaussian  distributed  velocity  component  to  drops. 


5.  Rainfall  Rate 

Rainfall  rate  determines  depth  of  accumulated  water  per  unit  time,  and  can  be 
derived  from  water  content  and  fall  speeds.  It  can  be  shown  that  the  a  cloud’s  water 
density  is 
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(2.38) 


M=^pw\D3N(D)dD 
®  o 

where  pw  is  the  water  density  of  a  drop. 

Turning  this  into  rainfall  rate  renders 

R  =  -]D3N(D)v(D)dD.  (2.39) 

6  0 


D.  RADAR  CROSS  SECTION  (RCS)  OF  RAIN 


1.  Cross  Section 

A  coherent  measure  of  radar  return  can  be  found  by  summing  the  electric  field 
contribution  from  each  point  scatterer  in  the  return  volume  of  interest.  Raindrops  can  be 
approximated  as  small  dielectric  spheres.  It  can  be  shown  (see  appendix  A)  that  the 
backscatter  field  of  a  small  sphere,  when  the  radius  is  much  smaller  than  a  wavelength 
( A  ),  is  given  by 


^  -  n  \v  I2  n6 

—  ^4  \KW\  D  , 


(2.40) 


where  D  is  the  drop  diameter,  K 


(£r~l) 
i£r+  2) 


and  £r  is  the  relative  permittivity. 


This  approximation  is  valid  as  long  as  the  wavelength  is  much  larger  than  the 
drop  circumference  and  the  drop  is  of  spherical  shape  and  is  referred  to  as  the  Rayleigh 
approximation  [3].  As  described  under  Physics  of  Precipitation,  drops  tend  to  flatten 
when  between  280  and  1000  pm  which  render  a  polarization  dependent  return.  This  can 
be  used  to  separate  different  sizes  of  precipitation  when  a  dual-polarization  radar  is  used 
to  obtain  the  ratio  of  two  orthogonally  polarized  returns  [3].  This  will  not  be  included  in 
the  initial  simulator  work. 

To  verify  the  return  predicted  by  the  simulator,  the  RCS  model  needs  to  be 
accurate.  Not  only  does  the  model  need  to  be  accurate  in  terms  of  magnitude,  the 
fluctuation  over  time  must  also  be  accurate  to  assure  the  model  will  produce  the  proper 
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variation  matching  the  theoretical  density  function.  For  the  magnitude  test  for  accuracy, 
Microwave  Studio  has  been  used  to  predict  the  scatter  from  spheres  with  properties  of 
raindrops.  Microwave  Studio  uses  the  Finite  Integration  Technique  to  calculate  RCS.  In 
Microwave  Studio  a  sphere  was  created  with  a  diameter  of  4  mm.  Conductivity  was  set 
to  0.0001  S/m  (pure  water)  [13].  The  set  up  for  all  simulations  was  in  most  aspects  the 
same  and  is  described  in  Figure  1  and  in  detail  in  appendix  B. 


z 


Figure  1.  Simulation  setup  using  Microwave  Studio,  single  drop.  kt  is 
the  incident  field  unit  vector  and  E  is  the  direction  of  the  electric  field. 

The  calculations  were  made  assuming  an  X-band  radar.  Although  drops  with 
D  >  4  mm  flatten  out  when  falling,  a  spherical  approximation  was  used  allowing  the 
Rayleigh  approximation  for  the  backscatter.  For  an  X-band  radar  with  frequency 

I  |2 

9.4  GHz  the  theoretical  backscatter  of  a  single  drop  with  \Km\  =  0.93  ,  and  a  diameter  of 
4  nun  is 

<7 b  =7r5^\Kmf  s(4xlQ  )  (9,4xl09)4  0.93m2 

b  ^4  I  ml  c4  {  ) 

<7b  =1.124x10 ~6m2  — >  -59. 5dBsm. 

The  result  using  Microwave  Studio  with  the  same  input  values  used  in  the  above 
is  shown  in  Figure  2  and  Figure  3. 
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Farfield  'farfield  (f=9.4)  [pw]'  RCS_AbsfTheta) 
90 


Figure  2.  Farfield  Scatter  from  a  sphere  of  4  mm  in  diameter  in  the  9 
plane  for  0  =  90°  and  0  =  270°. 


Farfield  'farfield  (f=9.4)  |pw]'  RCS_AbsfTheta) 


Figure  3.  Farfield  Scatter  from  a  sphere  of  4  mm  in  diameter  in  the  9 
plane  for  0  =  0  and  0  =  180° . 
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In  Figure  2  the  level  of  the  farfield  scatter  is  equal  in  all  directions,  which  is 
expected  considering  the  electric  field  direction  (x-direction)  and  the  view  of  observation. 
The  pattern  is  similar  to  a  dipole  where  the  farfield  scatter  has  its  maximum  in  the 
orthogonal  direction  from  the  electric  field  [15].  In  Figure  3  the  zero  directions  are  inline 
with  the  polarization  of  the  incident  filed,  which  occurs  when  0  =  0  and  0  =  180°  with 

e  =  90° . 


2.  Array  Factor 

Multiple  scatterers  can  be  modeled  as  an  array  of  scatterers  even  though  they  are 
not  identical.  One  must  realize  that  calculating  the  total  coherent  RCS  of  multiple 
scatterers  is  not  a  matter  of  summing  the  RCSs  since  RCS  is  acquired  by  summing  the 
scattered  electric  field  which  then  is  squared  to  get  the  RCS.  A  simple  approach  is  to 
calculate  the  electric  field  return  from  every  scatterer  and  then  add  the  contributions 
vectorally  [6].  The  phase  change  due  to  difference  in  location  using  a  far  field 
approximation  is  given  by 

e2Jks°" ,  (2.41) 

where  g0n  represents  the  path  difference  for  every  drop  with  reference  to  x0n ,  y0n ,  z0n  and 
is  given  by  (see  Figure  4) 

gon  =  x0n  sin  e  cos  </>  +  y0n  sin  e  sin  0  +  z0n  cos  0 .  (2.42) 
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X 

Figure  4.  Far-field  path  difference  due  to  shift  from  origin. 
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X 


Figure  5.  Simulation  setup  using  Microwave  Studio,  two  drops,  k.  is  the 
incident  field  unit  vector  and  E  is  the  direction  of  the  electric  field. 

The  value  is  calculated  only  for  a  direct  reflection  not  taking  into  account  the 
cos2  0  dependence  due  to  angle  of  incidence.  The  two  drop  setup  using  D  =  4  mm  and 
d  =  A  rendered  results  presented  in  Figure  6  and  Figure  7.  In  Figure  6  the  farfield  scatter 
pattern  show  a  maximum  in  the  orthogonal  direction  from  the  electric  field,  as  expected 
since  the  electric  field  was  along  the  y-axis  (0  =  90°  and  0  =  210°  with  0  =  90°).  In 
Figure  7  zero  directions  show  for  0  =  0  and  0  =  180°  with  0  =  30°  and  0  =  150° ,  which  is 
expected  since 


sin(Nkd  sin<?) 

= 

sin 

\ 

2  —  A  sin  0 
A  j 

sin  [4n  sin  0) 

sin  [kd  sin<9) 

sin 

f  2k  ,  .  3 

U  ) 

sin(2^sin6>) 

So  for  0  =  30°  and  <9  =  150°  the  expression  yields  zero.  Comparing  the  calculations, 
Equation  (2.45),  with  the  results  from  the  Microwave  Studio  simulations  reveals  a 
difference  of  only  0.1  dB,  which  can  be  the  effects  of  mutual  coupling.  The  max  values 
are  called  Bragg  lobes  and  their  locations  are  dependent  on  the  ratio  of  d  /  A ,  their 
location  with  respect  to  each  other,  and  the  measuring  point. 
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Faifield  ‘farfield  ((=9.4)  (pw)'  RCS_Abs (Theta) 


Figure  6.  Farfield  Scatter  from  two  spheres  of  4  mm  in  diameter  in  the 
0  -plane  for  0  =  90°  and  0  =  270°. 


Farfield  'farfield  (f=9.4)  (pw)'  RCS_Abs (Theta) 
90 


Figure  7.  Farfield  Scatter  from  two  spheres  of  4  mm  in  diameter  in  the 
0  -plane  for  0  =  0  and  0  =  180°. 
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The  question  of  mutual  coupling  was  addressed  using  Microwave  Studio.  The 
scattering  pattern  from  the  illumination  of  first  two  and  then  three  spheres  was  examined, 
this  time  for  different  separation  distances  (Figure  8).  The  values  plotted  are  the  results 
using  the  same  drop  setup  as  in  Figure  5,  i.e.  D-4  mm  ,  0  -  0  and  (j)  =  0  .  The  measures 
indicate  coupling  effects.  However,  the  software  is  not  optimized  to  calculate  low 
conductivity  targets,  as  raindrops  are.  When  the  space  between  the  drops  increases  the 
program  adds  noise  as  the  numerical  grid  becomes  significant  (see  Appendix  B). 
Although  mutual  coupling  is  likely  to  appear,  available  tools  could  not  provide  guidance 
how  to  provide  realistic  implementation. 


RCS  for  different  distances  between  drops 


Figure  8.  Monostatic  Back  Scatter  of  two  and  three  drops  of  4  mm 
diameter  when  spacing  is  varied  from  0.1  A  to  2  A. 

3.  Averaging  over  Time 

The  weather  radar  return  signal  is  composed  of  the  scattered  electric  field  from  all 
drops  in  the  resolution  cell  (see  Figure  9).  The  total  scattered  field  of  the  drops  in  the 
resolution  cell  is  given  by 
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(2.47) 


N 

E{r,0^)-a  =  Y,^eM'a  (2.48) 

tl= 1 

where  I//  is  given  by  i//n  =2 kg0n  and  is  a  function  of  r,  6  and  </>  (Equation  (2.43))  [6] 
and  a  is  a  unit  vector  giving  the  direction  of  the  scattered  field. 


z 


X 


Figure  9.  Radar  resolution  cell  used  in  Simulator.! 

It  can  be  argued  that  y/  is  a  random  variable  uniformly  distributed  on  the 
interval  [0,  2tt]  thus,  expected  RCS  will  be 


1  Azimuth  and  Elevation  angles  are  reversed  compared  to  the  mathematical  conventions. 
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(2.49) 


N 

(P total  )  —  CT„ 
n= 1 

where  (  )  denotes  the  expected  value.  So  to  get  the  average  power  estimate  a  number  of 
independent  samples  need  to  be  gathered  and  averaged. 

4.  Absorption  and  Attenuation 

Apart  from  the  fraction  of  the  incident  electric  field  that  rain  scatters,  rain  also 
absorbs  energy  which  causes  attenuation.  This  has  greater  effect  for  higher  frequencies, 
which  imposes  limits  on  the  possibilities  of  getting  higher  angular  resolution  by  using 
short-wavelength  radars.  The  effects  of  attenuation  can  be  used  in  bi-static  setups 
estimating  the  weather  parameters  based  on  the  level  of  attenuation.  Bi-static  setups  and 
effects  of  attenuation  will  not  be  taken  into  account  in  the  initial  simulator  work. 

E.  MATHEMATICAL  MODELS 

This  section  describes  the  basic  formulas  that  will  be  used  to  model  the 
parameters  required  in  the  simulator  to  generate  the  return  signal.  The  models  for  each 
parameter  are  defined  in  Table  1.  Although  choosing  Equation  (2.30)  as  ground  truth 
will  probably  create  too  high  estimates  of  reflectivity  [3],  this  will  be  the  initial  reference 
as  the  simulator  development  starts. 


26 


Table  1.  Mathematical  models  for  Weather  Estimators. 


Parameter 

Input  or  estimator 

Ground  Truth 

Drop  Size  Distribution, 

DSD 

N (D)  =  N0  exp(-AD), 

A  =  4.1 

N0  =8x1  O' m  mm  1 

Reflectivity  Factor,  Z 

1024i?2ln2/l2 

“  ~  Pfi.G'ClOdx3  |JC„f 

Z  =  ]D*N(D)dD  =  6-N° 

0  (4.  ITT0'21) 

Average  Radial  Velocity 

—  1  XT  \u  7/*  t 

V  =  iV-14* 

Compare  with  input  value 

Doppler  Spectrum  Width 

,  IM) 

O’2  =  M  ,  &  = 

v  At 

Compare  with  input  value 
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III.  SIMULATOR  DESIGN  CONSIDERATIONS  AND 
DEVELOPMENT  OF  SIMULATOR 


Having  described  the  required  parameter  models  in  the  previous  chapter  the 
design  and  development  of  the  simulator  will  commence.  As  stated,  the  simulator  should 
allow  for  inputs  of  rain  rate,  standard  radar  values  (e.g.  frequency,  azimuth  and  elevation 
angles,  and  pulsewidth),  wind  speeds  and  angles,  and  also  velocity  spread.  Output  signal 
parameter  estimators  of  interest  are  average  power,  average  Doppler  velocity  and 
Doppler  spread.  The  output  should  preferably  be  in  fonn  of  I  and  Q  return  voltage 
making  it  possible  to  post  process  the  data  to  obtain  the  estimate  of  average  power, 
Doppler  velocity  and  Doppler  spread. 


A.  SCOPE 

To  provide  a  tool  to  answer  the  research  questions  specified  in  the  introduction 
chapter,  the  following  is  required: 

1 .  A  physical  model  of  a  spatial  region  of  raindrops  must  be  created  and  later 
tested  and  verified  against  real  life  measures  or  existing  models, 

2.  A  model  for  realistic  motion  of  raindrops  that  are  weather  system 
dependent  must  be  created  where  at  least  an  interface  between  weather 
models  and  the  simulator  must  be  specified, 

3.  A  method  of  implementing  spatial  and  temporal  sampling  schemes  to 
study  the  output  for  constant  or  stepped  frequency  sampling, 

4.  A  method  to  study  the  use  of  pulse  compression  and  range  averaging  as  a 
means  of  rapidly  obtaining  independent  samples, 

5.  Displays  of  pedagogical  interest  to  facilitate  the  presentation  and  testing  of 
research  issues  mentioned, 

6.  A  study  of  the  perfonnance  of  estimators  for  the  first  three  Doppler 
moments, 

7.  A  study  of  the  utility  of  estimators  for  higher  order  Doppler  moments. 
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B.  LIMITATIONS 

The  initial  design  will  not  include: 

•  Dual  polarization 

•  Absorption/Attenuation 

•  Losses  and  noise 

•  Mutual  coupling  and  multiple  reflections 


C.  DEVELOPMENT 

Having  a  prototype  from  M.  Gosset  [1]  provided  valuable  ideas  and  experience 
developing  the  simulator.  The  method  used  to  develop  the  final  version  is  described  in 
Figure  10. 


Rewrite  Modules 
& 

Construct  Verification  Programs 


French  MATLAB  code 
Module  based  (functions) 


3  KJckjsj  jkl;j 

4  For  fnjkf 


7  JHffhjkhHKhJH 


Translation 
& 

Interpretation 


2  Gdjha  Gdhnmbd: 

3  KJckjsj  jkl;j 

7  JHffhjkhHKhJH 


Simulation 

Speed 

& 

GUI 


N(D)  =  N0  exp(-AD), 
A  =  4.1/?_0-2i/m/m“1, 

N0  =8xl03  m~3mm~' , 

R  =  ^DiN(D)v(D)dD 


Figure  10.  Plan  for  development  of  Simulator  using  the  French  version. 
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1.  Translation  and  Interpretation 

The  French  code  consists  of  a  main  program  that  calls  for  functions  that  performs 
different  tasks  to  support  the  main  program.  The  functions  are2: 

•  Size 

•  Approximation 

•  Marshall-Palmer 

•  Resolution  Volume 

•  Number  of  drops 

•  Initial  Position 

•  Summation 

•  Movement 

•  Atlas  Ulbrich 

The  input  parameters  are  shown  in  Table  2. 


Table  2.  Simulator  input  parameters 


Radar  Data 

Weather  Parameters 

Transmit  Power 

Drop  size  diameter  limits 

Max  Antenna  Gain 

Drop  Diameter  resolution 

Range  to  resolution  cell 

Rain  rate 

One-way  3  dB  Beamwidth 

Rain  angles 

Pulsewidth 

Wind,  speed  and  direction 

Frequency 

Pulse  Repetition  Frequency  (PRF) 

Elevation  and  Azimuth 

Number  of  pulses  for  integration 

2  The  code  was  translated  by  Professor  Monique  Fargues  of  NPS. 
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a.  Function  Description 

The  main  program  starts  by  calling  the  Size-function  that  generates  a 
discrete  Diameter  vector  based  on  the  diameter  interval  and  resolution  of  interest.  The 
function  calls  the  Approximation-function  that  calculates  the  RCS  for  each  diameter 
using  Equation  (2.40).  The  Marshall  Palmer-function  is  called  to  calculate  the  number  of 
drops  for  each  diameter  using  Equation  (2.35).  The  resolution  is  multiplied  with  the 
number  of  drops  per  diameter  and  stored  in  a  matrix  together  with  the  Diameters  and 
RCSs.  The  Size-function  also  has  a  part  that  uses  a  pre-set  threshold  to  take  away  the 
smaller  drops  that  account  for  10  %  of  the  total  RCS.  The  matrix  created  in  the  Size- 
function  consists  of  Diameter,  RCS,  Number  of  drops,  total  RCS  per  diameter,  and 
Reflectivity  Factor  structured  as 


'  A 

D„  ' 

v(a) 

...  N(D„) 

A-tf(A) 

l  Zi 

ZN  j 

where  the  first  values  are  set  by  the  threshold  and  N  by  the  upper  Diameter  limit. 


(3.1) 


The  main  program  then  calls  the  Resolution  Volume-function  to  generate 
the  resolution  volume  using  the  radar  data.  The  function  adds  a  margin  to  the  resolution 
volume  so  that  the  drops  can  pass  through  the  resolution  volume  as  time  passes  between 
pulses.  Both  the  borders  for  the  increased  resolution  volume,  henceforth  called  the  Box, 
and  the  resolution  volume  are  defined  for  later  use.  The  main  program  uses  the  data  to 
feed  the  Number  of  drops-function  that  calculates  the  number  of  drops  considering  the 
volume  of  the  Box  and  changing  the  values  in  the  Matrix. 


Using  the  Initial  Position- function,  the  simulator  uniformly  distributes 
every  drop  inside  the  Box  using  the  data  from  the  Matrix.  The  positions  are  stored  in 
files  unique  for  each  diameter  using  spherical  coordinates  to  simplify  comparison  when 
determining  whether  a  drop  is  inside  the  Resolution  Volume  or  not. 


The  main  program  then  starts  the  simulation  for  the  number  of  pulses 
specified.  For  all  iterations  the  Summation-function  calculates  the  return  voltage  and  the 
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return  power  taking  into  account  the  position  and  RCS  of  each  drop.  The  radar  cross 
section  is  calculated  summing  the  scattered  electric  field  from  each  rain  drop  using 

N, 

E'0lai  {h )  =  ZZ  yl°D.  •  ej2kr‘[,k) 

*  (3.2) 

where  tk  is  the  sample  time  and  s  is  the  Diameter  vector. 

The  results  are  stored  in  different  files  depending  on  what  data  it  is. 
Before  next  iteration  the  main  program  calls  the  Movement-Function  to  move  the  drops 
according  to  wind  and  fall  velocities.  To  get  the  estimated  tenninal  velocities  for  each 
diameter  the  Atlas  Ulbrich-function  [12]  is  used,  which  is  based  on  Equation  (2.37). 

The  output  of  the  simulator  consists  of  four  plots  that  show  the  Box  with 
all  initial  drops  plotted,  the  Resolution  Volume  with  all  positions  the  drops  have  during 
the  simulation  plotted,  power  return  levels  for  every  pulse,  and  the  phase  of  every  drop 
and  every  pulse. 

b.  Verification 

To  evaluate  the  accuracy  of  the  simulator,  verification  programs  were 
created  to  estimate  errors  of  the  modules  of  the  simulator.  The  modules  to  examine  were 

•  Number  of  drops 

•  Average  power  estimation 

•  Doppler  estimation 

•  Doppler  spread  estimation 

(1)  Number  of  drops.  The  original  French  program  used  a 
diameter  interval  between  0.05  mm  and  5  mm  with  a  resolution  varying  from  0.08  mm 
for  the  smaller  drops  (up  to  1  mm)  to  0.2  mm.  Using  the  original  rain  rate  set  to  50  mm 
renders  a  histogram  as  seen  in  Figure  1 1 .  The  threshold,  that  removes  the  smaller  drops 
that  account  for  10  %  of  the  total  RCS,  is  also  indicated  in  the  same  figure. 
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Drop  Size  [mm] 


Figure  11.  Drop  Size  Distribution  comparison  for  French  simulator.  The 
volume  is  set  to  1  m3  and  the  rain  rate  to  50  mm. 

Eliminating  the  smallest  drops  will  effect  the  fluctuation  of  the 
return  signal  as  well  as  the  peak  values,  which  in  turn  affects  the  average  value.  This  is 
an  effect  of  the  array  factor  previously  described  that  can  make  several  small  drops 
generate  high  return  values  when  in  phase.  Also  the  choice  of  upper  limits  will  affect  the 
level  of  power  return  especially  for  higher  rain  rates,  which  will  be  examined  in  a  later 
section.  These  values  are  plotted  without  making  the  number  of  drops  to  integer  values. 
The  simulator  uses  CEIL  to  round  off  the  number  of  drops  making  the  integer  equal  to  or 
higher  than  the  original  value.  Simulations  show  that  the  effect  becomes  significant  for 
rain  rates  below  40  mm/hr,  creating  too  high  return  power  values  (Figure  12,  RCS  using 
original  French  simulator).  To  produce  values  closer  to  what  the  integral  form  generates 
using  Equation  (2.30),  the  span  of  diameters  needs  to  be  widened  and  the  round  off  tool 
exchanged. 

(2)  Average  Power  Estimation.  To  verify  the  average  power 
return,  a  simplified  simulator  was  built  that  didn’t  include  the  movement  of  the  drops. 
The  movement  was  simulated  by  randomly  re-plotting  the  drops  for  all  iterations.  Using 
the  knowledge  about  the  effects  of  choosing  limits  for  the  drops  size,  a  simulation  was 
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performed  verifying  the  levels  for  average  return  power.  Also  the  statistical  distribution 
of  the  return  signal  was  examined  to  assure  it  agreed  with  basic  theory  [1],  [2],  The 
simulator  was  run  for  1000  samples  and  for  rain  rates  between  10  and  120  mm/h  and  a 
wavelength  A  =  0.1  m. 

The  results  (Figure  12)  show  that  using  the  original  French 
simulator  makes  the  relation  between  RCS  and  rain  rate  almost  linear  while  the  new 
generator  follows  the  integral  form  (Equation  (2.30)).  The  new  generator  uses  a  wider 
span  of  diameters  (up  to  15  mm),  does  not  apply  a  threshold  and  uses  ROUND  instead  of 
CEIL  making  the  value  closer  to  the  real  value. 

RCS  density  of  1  m3  of  drops. 
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Figure  12.  Average  RCS  density  comparison  using  different  drop 
generators.  1000  samples  were  used  and  compared  with  the  result  using  the  integral 
form  of  reflectivity  to  calculate  RCS  density.  A  =  0.1  m. 


The  new  generator  still  produces  values  with  errors  up  to  2  dB 
even  with  sample  sizes  as  large  as  1000  samples.  To  examine  the  cause  of  this  error 
another  version  of  a  drop  generator  was  constructed  taking  into  account  the  error  caused 
by  the  rounding  off.  The  generator  keeps  track  of  the  error  and  adds  it  at  the  end  of  the 

simulation.  This  also  provided  an  idea  of  how  to  be  able  to  improve  the  simulator  speed 
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by  reducing  the  number  of  drops  while  still  retrieving  accurate  RCS  and  statistical 
behavior.  Reducing  the  number  of  drops  by  dividing  by  a  constant  and  after  the 
simulations  multiply  it  back  into  the  calculations  to  retrieve  reflectivity,  number  of  drops, 
and  rain  rate  is  valid  since  all  these  operations  are  performed  by  integrating  over  the 
diameter.  Treating  the  error  like  a  constant  results  in  a  max  error  of  about  0.4  dB,  an 
error  that  seems  to  increase  as  the  rain  rate  increases  (see  Figure  13).  Since  the  simulator 
replaces  the  integration  by  a  summation,  errors  will  arise.  Using  Equation  (2.30)  renders 
a  Z-value  of  54.1  dBZ,  for  a  rain  rate  of  100  mm/h.  When  using  the  same  form  of 
summation  the  simulator  applies,  the  same  estimation  is  54.5  dBZ,  which  explains  the 
error  of  0.4  dB. 
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Figure  13.  RCS  density  comparison  using  new  drop  generator  with  error 
correction  and  reduction  factor.  A  -  0. 1  m. 

To  get  the  average  value  of  the  return  power,  several  independent 
samples  are  used.  The  RCS  of  a  raindrop  was  derived  using  the  radiation  integral,  which 
can  be  used  when  the  return  complex  voltage  is  calculated.  Since  I  and  Q  are  random 
variables,  independent  and  none  of  them  dominant,  the  central  limit  theorem  states  that 
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their  sum  tends  to  be  Gaussian,  as  long  as  the  numbers  of  samples  are  large  enough.  The 
amplitude  of  the  return  voltage  i.e.  the  absolute  value  of  the  I  and  Q  voltage  return, 

\V\  =  (l2  +  Q2)Y  (3.3) 

is  therefore  expected  to  have  a  Rayleigh  distribution  [3],  Running  the  drop  generator  for 
10,000  samples  renders  results  shown  in  Figure  14  and  Figure  15. 


RCS  density  of  Irrr  of  drops.  Rainrate=  50.  Number  of  drops=  4454.  0=  0,  0=  1.5708 
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Figure  14.  Plot  of  10,000-sample-run  of  drop  generator. 
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Figure  15.  Rayleigh  fit  of  data  from  10,000-sample  run  of  drop  generator. 


In  Figure  14  the  average  value  of  RCS  is  lower  than  the  result 


using  the  integral  form.  This  was  explained  earlier  when  the  results  from  summing  over 
all  diameters  was  compared  with  the  integration  from  zero  to  infinity.  The  Rayleigh  fit  in 
Figure  15  follows  the  expected  distribution.  The  results  were  retrieved  by  plotting  a 
histogram  of  the  10,000-sample-run.  Both  the  histogram  and  the  tit  were  normalized  to 
allow  for  comparison. 


(3)  Doppler  estimation.  The  module  that  moves  the  drops  takes 


into  account  the  terminal  velocities,  determined  by  Equation  (2.37)  with  an  associated  fall 
angle,  and  wind  speeds  in  x,  y,  and  z-dimension.  The  spread  of  terminal  velocities 
together  with  the  fall  angle  will  contribute  to  the  spread  of  the  Doppler  spectrum.  To 
evaluate  the  mean  Doppler  frequency  a  pulse  pair  algorithm  is  used  to  record  the  phase 
change  over  time  using  Equation  (2.32)  and  comparing  the  phase  extracted  from  two 
consecutive  returns  pulses.  Applying  this  on  the  drop  generator  and  the  movement 
module  enables  an  estimate  of  the  mean  frequency.  A  simulation  was  run  for  30  pulses 
(29  pulse  pairs)  with  a  PRF  of  2000  Hz.  Wind  speed  was  set  to  zero  and  the  rain  fall 
angle  to  30°  with  respect  to  the  z-axis  (for  detailed  specifications  see  appendix  C).  The 
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simulation  showed  an  average  Doppler  frequency  fd  =-57.3  Hz  or  vr  =-2.85  m/s  in  x- 

direction  (radial  direction)  and  standard  deviation  of  o f  =  26.8  Hz,  which  represents 

<JV  =1.34  m/s  (see  Figure  16).  The  spread  was  induced  by  the  differences  in  terminal 

velocities  for  different  drops  sizes.  To  estimate  the  accuracy,  calculations  were  made 
using  a  2  mm  drop  for  which 

v(D)  ~  38 6.6D°'67  =  386.6 ■  0.002° 67  =  7.89[/» / s\ 

Taking  the  fall  angle  into  account  renders 

vx_dir  =v-sin(30o)cos(0)  =  -3.0l[m/,y] . 

To  further  verify  that  the  simulator  generates  accurate  results,  the 
fall  angle  was  removed  and  only  wind  in  the  x-direction  was  used.  The  same  inputs  were 

used  adding  wind  of  -5  m/s  rendering  fd  =-105.1  Hz,  which  corresponds  to 
vr  =-5.3  m/s  as  seen  in  Figure  17.  The  Doppler  spread  was  6 f  =7.9  Hz  ( <JV  =  0.4 

m/s).  Note  that  the  first  spectrum  (Figure  16)  is  skewed  towards  zero  indicating  that  the 
smaller  drops,  constituting  the  majority  of  the  drops,  bias  the  result.  Also  the  limited 
number  of  samples  (29  pulse  pairs)  can  have  effects  on  the  shape  of  the  histogram. 
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Figure  16.  Doppler  Frequency  Histogram  applying  terminal  velocities  and 
30°  fall  angle.  30  pulses  and  29  pulse  pairs. 
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Figure  17.  Doppler  Frequency  Histogram  applying  -5  m/s  wind  in  -x- 
direction.  30  pulses  and  29  pulse  pairs. 
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(4)  Doppler  spread.  The  simulation  run  above  showed  a  standard 
deviation  of  1.18  m/s,  which  is  a  measure  of  how  much  the  average  radial  velocity  of  all 
scatterers  of  each  pulse  return  spread  over  all  samples.  The  data  was  fitted  to  a  normal 
density  function,  which  is  expected  for  a  power  spectrum  [3],  to  explore  if  the  data  could 
be  considered  Gaussian.  The  second  run  fits  better  with  the  Gaussian  plot,  which  can  be 
explained  by  the  randomness  of  the  process  different  from  the  first  simulation  where  the 
spread  is  distorted  by  the  differences  in  drop  sizes  and  therefore  terminal  velocities.  To 
further  examine  the  power  spectrum  new  algorithms  were  implemented  as  described  in 
the  upcoming  section. 

2.  Adding  New  Modules 

The  verification  programs  have  been  used  to  validate  the  simulator  but  there  are 
still  improvements  that  can  be  made  to  achieve  better  accuracy  and  realism.  First  to 
recapitulate  what  has  been  accomplished  up  until  this  point: 

A  drop  generator  that  provides  the  number  of  drops  given  by  the 
distribution  of  interest  based  on  rain  rate. 

An  algorithm  that  calculates  the  complex  scattered  voltage  return  based 
on  the  position  of  each  drop  within  the  radar  resolution  volume. 

A  module  that  introduces  movement  to  the  drops  concerning  fall  speed 
with  an  associated  angle,  and  wind  speed. 

An  algorithm  that  calculates  the  doppler  frequency  of  the  moving  drops 
based  on  phase  change  between  pulse  pairs. 

An  error  correcting  algorithm  that  also  can  be  used  to  reduce  the  number 
of  drops  used  in  a  simulation  to  make  it  faster. 


a.  Range  Weighting 

To  provide  a  more  realistic  level  of  the  return  signal  value,  a  range 
weighting  function,  based  on  the  radar  range  equation,  was  built  and  implemented.  As 
mentioned  earlier,  effects  of  absorption  and  other  losses  will  not  be  considered,  leaving 
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only  the  effects  of  difference  in  range.  To  account  for  range  effects,  Equations  (2.21), 
(2.24),  and  (2.48)  were  used  rendering 
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where  a  is  a  unit  vector  giving  the  direction  of  the  scattered  field.  For  the  total 
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where  y/n  =  2kg  0n  and  g0n  is  given  by  Equation  (2.42). 


b.  Correlation 

Another  estimate  of  the  doppler  from  frequency  the  return  electric  field 
can  be  obtained  by  taking  the  discrete  Fourier  transform  (here  FFT)  of  the  signal  samples 
using  the  complex  voltage  return.  The  Doppler  Power  Spectrum  is  given  by 

s(/)=|z(/)|T,/m 

m- i  (3-6) 

m= 0 

where  M  is  the  number  of  samples.  An  estimate  can  be  retrieved  by  taking  the  Fourier 
transfonn  of  the  autocorrelation  function 

*(/)  =  —  Z  V*(m)V(m  +  l).  (3.7) 

M  m-Q 
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Computing  the  autocorrelation  function  and  the  FFT  for  the  output  in  the 

example  on  pg.  39  produces  the  same  average  frequency  estimate,  /  = -105.1  Hz  as 
before  (see  Figure  18). 
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Figure  18.  Doppler  spectrum  using  FFT  and  the  Fourier  transform  of  the 
autocorrelation  function.  (They  overlap  each  other).  Peak  gives  average  frequency 

estimate  /  =-105.1  Hz. 

Both  autocorrelation  and  FFT  preserve  amplitude  information,  which  has 
consequences  interpreting  average  Doppler  using  these  functions.  Larger  drops  will 
affect  the  power  spectrum  to  a  greater  extent  than  smaller  since  drop  RCS  is  proportional 
to  Db .  Thus,  any  difference  in  the  radial  velocity  of  large  and  small  drops  will  result  in 
the  power  spectrum  being  skewed  towards  the  Doppler  associated  with  the  larger  drops. 
Using  the  pulse  pair  algorithm  (Equation  (2.13))  removes  the  amplitude  information 
producing  a  more  “true”  average  doppler  frequency  or  radial  velocity. 

c.  Frequency  Stepping 

To  accomplish  rapid  volumetric  scanning  there  are  several  methods 
proposed.  One  method  to  achieve  de-correlation  is  by  shifting  in  frequency.  A  fast 
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scanning  scheme  to  speed  up  the  acquisition  of  independent  samples  was  proposed  by 
Doviak  and  Zrnic  [3].  They  suggested  a  method  where  a  pulse  pair  is  sent  with  one 
frequency  assuring  dependent  doppler  measures,  then  the  frequency  is  changed  a 
minimum  of  1/r,  for  de-correlation  to  send  the  next  pulse  pair.  Implementing  this 
scheme  into  the  simulator  requires  modules  stepping  between  more  than  one  PRF  also 
taking  into  account  the  differences  in  wavelengths  when  computing  the  meteorological 
signal  parameter  estimates.  In  Figure  19,  a  three  frequency  stepping  scheme  is 
introduced  that  transmits  a  pulse  pair  at  frequency  1,  and  listens  for  the  time  set  by  the 
Doppler  Pulse  Repetition  Time,  PRTi,  before  changing  to  frequency  2.  Before  the 
second  frequency  can  be  transmitted  the  frequency  oscillator  needs  a  settling  time,  here 
called  delay.  After  all  three  frequencies  have  been  transmitted  the  radar  returns  to 
frequency  1 .  The  waiting  time  will  be  dependent  on  the  de-correlation  time  for  frequency 
1,  which  in  this  case  is  PRT2  rendering  a  true  time  delay  before  resending  frequency  1  of 

t delay  =  PRT2  ~  N fy  (PRTX  +  Delay) 

where  Nf  equals  the  number  of  frequencies  used. 

Pulse  1  Pulse  2 

Frequency  1 


Frequency  2 


Frequency  3 


Figure  19.  Frequency  scheme  using  stepped  frequencies  to  obtain  de- 

correlation. 


44 


d.  Pulse  Compression 

A  benefit  of  Pulse  Compression  lies  in  the  capability  of  averaging  over 
wider  range  where  the  gain  in  range  resolution  provides  independent  samples  that  can  be 
used  to  estimate  the  average  power  for  the  volume  of  the  wider  pulse.  Implementing 
chirped  pulse  compression  in  the  simulator  calls  for  the  ability  to  step  in  frequency  as 
well  as  range.  A  scheme  was  developed  enabling  stepped  frequency  and  range  (see 
Figure  20),  where  the  Resolution  Cell  is  divided  into  smaller  range  bins  defined  by  the 
sub  pulse  width  and  the  number  of  frequency  steps. 
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Figure  20.  Pulse  Compression  implementation. 


When  the  frequency  has  been  stepped  through  the  whole  pulse,  the  drops 
are  moved  in  accordance  to  the  chosen  PRT.  The  process  then  starts  the  frequency  and 
range  stepping  again,  collecting  returns  from  the  next  pulse.  The  total  power  from  each 
pulse  is  summed  over  all  frequencies  and  then  averaged  for  the  number  of  range  bins 
used. 


e.  Plot  of  Phasors  for  Phase  Change  Comparison 

To  visually  illustrate  the  phase  change  for  pulse  pairs  over  time,  30  pulse 

pairs  are  presented  in  polar  plots  as  phasors  (see  example  Figure  21).  The  current  version 

of  the  simulator  can  only  plot  30  pulse  pairs  but  uses  all  simulation  pulse  pairs  to 
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calculate  the  average  radial  velocity  estimate.  The  estimated  radial  velocities  retrieved 
from  each  pulse  pair  are  plotted  in  a  Histogram  with  a  Gaussian  fit  for  comparison,  as 
presented  in  Figure  22.  A  three  frequency  stepped  scheme  simulation  (pg.  44)  produced 
an  average  radial  velocity  estimate  of  vr=-7.38  m/s  and  a  standard  deviation  of 
<7V  =  4.68  m/s  for  a  single  run  of  51  pulses  calculated  using  Equations  (2.32),  (2.33)  and 
(2.34)  (see  Figure  22). 


Radial  Velocity  Histogram  with  Gaussian  Fit 


Radial  velocity  [m/s] 


Figure  22,  Radial  Velocity  Histogram  with  Gaussian  Fit  for  51  pulse 
pairs.  Input  velocity  was  -10  m/s  with  a  velocity  spread  of  4  m/s. 
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f.  Adding  Spread  to  the  Doppler  Spectrum 

As  described  in  a  previous  section,  a  spread  of  the  velocity  spectrum  needs 
to  be  implemented  to  make  differences  in  correlation  possible.  Applying  a  Gaussian 
random  component  to  the  wind  velocity  to  introduce  a  controllable  value  of  velocity 
spread  was  tested.  This  enables  verification  of  the  theory  of  the  spread  contributing  to 
de-correlation.  Results  from  running  simulations,  again  using  a  PRF  of  2000  Hz  and  no 
fall  angle,  for  different  spreads  applied  can  be  observed  in  Figure  23.  The  spread  applied 
was,  from  left  to  right,  <Jv  =  4  m/s,  <Jy  =  3  m/s,  and  er,  =  2  m/s,  and  the  average 

simulator  velocity  spread  estimates  were  ov  =  4.05  m/s,  <7V  =  2.94  m/s,  and 
<7v  =  2.00  m/s  again  calculated  using  Equations  (2.32),  (2.33)  and  (2.34). 


Doppler  Frequency  [Hz]  Doppler  Frequency  [Hz]  Doppler  Frequency  [Hz] 

Figure  23.  Test  of  Doppler  spread,  (a)  input  spread  was  4  m/s,  output 
4.05  m/s.  (b)  input  spread  was  3  m/s,  output  2.94  m/s.  (c)  input  spread  was  2  m/s, 

output  2.00  m/s. 

The  best  de-correlation  results  were  accomplished  when  the  velocity 
spread  was  applied  in  two  steps;  first,  a  spread  representing  the  spread  of  wind  velocity 
over  time,  second,  a  spread  representing  the  spread  due  to  differences  in  drop  sizes  and 
other  effects  that  will  make  the  velocity  of  each  drop  vary  over  time.  The  fact  that  drops 
fall  out  of  the  resolution  volume  and  new  drops  join  adds  an  uncontrollable  spread 
mechanism  that  can  be  examined  by  running  several  consecutive  runs  with  the  same 
input  data.  As  the  number  of  samples  increases,  the  spread  of  the  average  velocity  is 
expected  to  become  smaller  [3],  which  can  be  used  for  comparison. 
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D.  SIMULATOR  DESCRIPTION 

The  simulator  was  built  in  modules  developed  to  enable  desired  functionalities. 
Initially  the  modules  were  separate  functions  that  a  main  program  called  for  execution, 
but  MATLAB  works  faster  when  only  one  long  script  is  used  rather  than  functions.  The 
inputs  are  separated  into  Radar  data  and  Rain  parameters  and  are  specified  in  previous 
sections.  Part  I  of  the  simulator  uses  input  Radar  data  to  generate  a  Beam  Resolution  cell 
approximation  and  a  Cube  with  an  added  margin  to  the  Beam  Resolution  cell  to  allow  for 
the  rain  drops  to  fall  the  through  beam  over  time.  Part  II  uses  the  Rain  Parameters  to 
generate  a  fundamental  T  Matrix  consisting  of: 

•  Diameter  Vector  and  Number  of  drops  in  each  diameter  interval  using  the 
Marshall-Palmer  exponential  approximation  for  drop  size  distribution. 

•  The  backscatter  RCSs  for  each  Diameter  size 

•  The  total  number  of  drops  of  each  diameter  per  m 

When  the  frequency  stepped  version  is  used  the  RCSs  are  placed  in  a  separate 
matrix,  since  RCS  is  frequency  dependent.  In  Part  III  all  drops  a  placed  randomly  in  the 
Cube  developed  in  Part  I  and  in  Part  IV  the  coherent  Electric  field  return  is  calculated  for 
every  pulse.  Between  pulses  all  drops  are  moved  in  accordance  with  their  velocities 
where  velocity  has  a  random  component  determined  by  <Jv .  In  Part  V  the  final 
calculations  are  completed  and  the  results  are  plotted. 


1.  Part  I 

To  create  the  radar  resolution  cell  based  on  the  Radar  data  the  volume  of  interest 
is  retrieved  by  calculating  the  volume  of  the  cone  segment  that  represents  the  resolution 
cell  (see  Figure  9).  Assuming  an  elliptical  beam,  the  resolution  volume  can  be 
approximated  as  [5] 


V_ 


n 


{  CT^ 


V  Z  J 


(3.8) 


To  allow  for  the  water  drops  to  fall  through  the  resolution  cell  over  all  pulses,  a 
margin  is  added  to  the  resolution  cell  creating  a  box  where  all  drops  are  initially  placed. 
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Both  volumes  are  defined  in  spherical  coordinates  to  allow  for  placement,  plots,  and  later 
comparison.  Also  the  volumes  are  registered. 


2.  Part  II 

This  part  produces  a  matrix  of  all  drops  that  contribute  to  the  total  reflectivity  of 
the  resolution  cell.  Using  the  diameter  limits  defined  in  Rain  parameters,  a  drop  size 
vector  is  created.  The  limits  can  be  chosen  to  create  results  close  to  the  integral  form  by 
pushing  the  upper  limit  of  diameter  to  15  mm  (see  results  Figure  12),  or  to  resemble  more 
realistic  values  setting  the  max  diameter  to  around  6  mm.  Depending  on  what  limits  are 
chosen,  the  actual  Z  must  be  computed  using  the  appropriate  model,  i.e.  chose  between 
Equation  (2.30)  and  (2.31).  The  simulator  allows  for  differences  in  resolution  for 
different  diameter  ranges  depending  on  what  one  wants  to  examine.  The  diameter  vector 
is  then  used  to  create  an  RCS  vector  applying  Equation  (2.40),  and  a  Drop  Size 
Distribution  (DSD)  vector  applying  Equation  (2.35)  (Marshall-Palmer).  To  improve  the 
speed  of  the  simulator  the  number  of  drops  is  scaled  down  by  a  factor  (Norm),  which  is 
used  later  to  rescale  the  output  values.  The  number  of  drops  is  rounded  off  to  nearest 
integer  using  the  MATLAB  function  ROUND.  Finally  a  matrix  is  created  to  store  all 
values  computed  in  Part  II 


f 

matrix  T  round  = 
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D, 


(3.9) 


WA)  •••  n(dn)) 

and  errors  due  to  rounding  off  are  calculated.  The  errors  accounted  for  are  RCS,  Z, 
Number  of  Drops,  and  Rain  Rate.  In  the  stepped  frequency  and  pulse  compression 
versions  the  RCSs  are  stored  in  a  separate  matrix 
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3.  Part  III 

The  drops  created  in  Part  II  are  placed  inside  the  Box  created  in  Part  I.  For  each 

diameter  every  drop  is  randomly  placed  inside  the  Box  using  a  uniform  distribution 
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function  and  its  position  is  stored  in  spherical  coordinates  in  a  position  matrix  saved  in  a 
.mat-  file 


f 


matrixpos  = 


\anglex 


where  angle  is  the  angle  off  boresight. 
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(3-11) 


4.  Part  IV 

In  Part  IV  the  Coherent  Electric  Field  return  is  calculated  for  every  pulse 
specified  in  the  Radar  Data.  For  every  diameter,  drops  that  are  positioned  inside  the 
Resolution  Cell  are  selected  and  based  on  matrices  (3.9)  and  (3.11)  the  weighted  electric 
field  return  is  calculated  using  Equation  (3.5)  and  summed  for  all  drop  diameters.  Also 
the  rain  rate,  for  the  specific  Resolution  Cell  is  calculated  using  Equation  (2.39) 
approximated  by  a  summation  over  all  drop  sizes.  After  the  calculations,  all  drops  are 
moved  to  a  new  position  by  applying  the  wind  vector  and  terminal  fall  velocities 
specified  in  Rain  parameters  using  Equation  (2.37).  All  drops  are  assumed  to  move  with 
the  velocity  of  the  wind  vector  with  a  Gaussian  spread  applied.  The  spread  can  be 
applied  equally  to  all  drops,  so  all  drops  move  with  the  same  velocity  between  pulses  and 
the  spread  is  applied  over  all  pulses,  or  applied  per  diameter  size  so  equally  sized  drops 
move  at  the  same  velocity  between  two  specific  pulses  and  the  spread  is  applied  over 
sizes.  A  third  possibility  is  to  apply  the  velocity  spread  on  all  drops  individually  creating 
a  true  spread  although  not  fully  realistic.  The  new  drop  positions  are  stored  in  matrixpos. 

In  the  stepped  frequency  version,  two  different  PRTs  are  used  together  with  a 
delay  described  in  previous  sections.  After  the  first  pulse  has  been  registered  in  terms  of 
I  and  Q  return,  using  frequency  1,  the  drops  are  moved  according  to  the  Doppler  PRT 
( PRT1 ).  The  phase  is  then  registered  completing  the  pulse  pair  before  moving  the  drops 
again,  this  time  using  the  Doppler  PRT  and  the  delay  marking  the  change  of  frequency. 
Frequency  2  and  frequency  3  are  then  applied  in  the  same  way  as  frequency  1.  When 
changing  back  to  frequency  1  the  simulator  uses  the  second  PRT  ( PRT _2 )  and  subtracts 
the  time  elapsed  due  to  the  other  frequencies  and  their  transmitted  pulse  pairs.  The  return 
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values  are  stored  in  a  matrix  making  sure  the  errors  and  the  final  estimates  account  for  the 
differences  in  wavelengths  and  range  between  the  returns. 

The  pulse  compression  version  steps  through  the  number  of  stretched  pulses 
specified  in  the  Radar  data.  For  every  stretched  pulse  the  chirp  is  approximated  by 
discrete  frequency  steps  using  the  frequency  step  size  and  the  number  of  steps  defined. 
Every  frequency  is  associated  with  a  specific  range  bin  defined  by  the  Resolution  Cell 
depth  and  the  number  of  frequency  steps.  After  the  return  values  for  all  frequencies  and 
range  bins  have  been  collected,  all  drops  are  moved  and  a  new  pulse  initiated.  The  time 
between  pulses  is  defined  by  the  PRT. 

5.  Part  V 

The  final  part  of  the  simulation  produces  the  weather  signal  parameter  estimates 
and  plots  the  results.  The  estimated  parameters  are 

•  Average  Power  estimate,  P,  retrieved  using  the  weighted  electric  field 
returns. 

•  Reflectivity  estimate,  Z ,  retrieved  from  the  average  power  estimate  using 
(2.28)  and  compared  with  Z  based  on  the  simulator  inputs  using  (2.27), 
evaluated  by  a  summation  and  Equation  (2.30),  evaluated  using  the 
parameters  N0  and  A  from  the  Marshall-Palmer  DSD. 

•  Average  doppler  frequency  estimate,  fd ,  retrieved  using  a  pulse  pair 
algorithm  that  extracts  the  phase-change  over  time  (see  Equation  (2.32)). 

•  Doppler  spectrum  width  estimate,  & f ,  retrieved  taking  the  standard  deviation 
of  the  doppler  frequency  data  (see  Equation  (2.33)  and  (2.34)). 

Here  the  errors  and  the  reduction  factor  are  put  back  in  the  calculations.  As  a 
comparison,  the  ground  truth  for  the  respective  parameter  is  derived  using  the  integral 
form  for  calculation  and  plotted  in  the  same  graphs  for  clarity.  Examples  of  output  plots 
are  presented  in  Figure  24  where  a  simulation  using  three  frequency  steps  was  evaluated. 
The  input  parameters  include  frequency  steps  from  3  GHz  to  3.002  GHz  with  a  1  MHz 
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frequency  step,  a  Doppler  PRF  of  5,000  Hz,  and  the  second  PRF  of  300  Hz.  Further,  the 
radial  velocity  was  -10  m/s  with  a  velocity  spread  of  4  m/s,  which  should  produce  a 


p  ~  0.246  ,  using  T  = - 

s  300 


s. 


The  average  Doppler  frequency  estimate  fd  =  -1 50  Hz 


(Figure  24  (a))  is  about  25  %  from  the  expected  Doppler  frequency  of  fd  =  -200  Hz. 

Only  5 1  sample  pairs  were  used  to  estimate  the  average  Doppler,  which  could  explain  the 
difference  between  the  estimate  and  the  true  value.  The  simulation  is  a  point  estimate, 
which  means  that  running  consecutive  trials  would  better  show  the  behavior  of  the 
velocity  estimate.  The  offset  can  also  be  the  effects  of  adding  new  drops  and  loosing 
others  in  the  resolution  cell  as  explained  in  previous  sections.  The  plot  of  Estimation  of 
Z,  (c)  presents  results  close  to  the  outcome  expected  using  the  mathematical  formulas  for 
comparison.  The  correlation  value  of  0.246  can  explain  why  the  value  of  the  estimate  Z 
is  so  close  to  the  integral  fonn,  as  de-correlated  samples  will  yield  an  estimate  with  a 
smaller  variance.  The  example  output  plot  of  correlation  coefficient,  (d),  although  close 
to  the  expected  value  from  the  above  described  simulation,  is  from  another  simulation 
running  only  one  PRF.  The  input  parameters  were: 


•  Wavelength,  X  0.1  m 


•  PRF  300  Hz 


•  Wind  -10  m/s  in  x-direction  (=radial  velocity) 

•  Velocity  spread  4  m/s 

The  expected  correlation,  using  Equation  (2.15),  is  p  ~  0.246  ,  which  is  very  close 
to  the  simulator  output  value. 
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(a) 


(b) 


Doppler  frequency  for  all  pulse  pairs  for  PRF=  5000  Hz 


Power  Return  for  all  51  pulses.  Frequency  stepped  from  3  GHz  to  3.002  GHz  with  1  MHz  steps. 


(C) 


Estimated  Z  per  pulse 


(d) 


Return  Voltage  Correlation 


Figure  24.  Example  output  plots  from  simulator,  (a)  Doppler,  (b)  power, 
(c)  reflectivity,  and  (d)  correlation  coefficient.  Note  that  (d)  is  from  a  different 

simulation  than  (a)  -  (c). 
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IV.  ANALYSIS  OF  APPLICATIONS  AND  DISCUSSION  OF 

RESULTS 


A.  IS  IT  POSSIBLE  TO  BUILD  A  WEATHER  RADAR  SIGNAL 

SIMULATOR  BASED  ON  A  PHYSICAL  MODEL  OF  A  SPATIAL 

REGION  CONTAINING  RAINDROPS? 

The  physical  model  built  is  based  on  the  Rayleigh  scattering  approximation  and 
the  number  of  scatterers  and  their  sizes  are  specified  using  the  Marshal-Palmer  drop  size 
distribution.  The  drops  are  approximated  by  spheres  and  are  unifonnly  distributed  in  a 
volume  defined  by  the  radar  data  and  the  return  signal  is  calculated  by  summing  the 
electric  field.  The  model  accounts  for  varying  terminal  velocities  for  different  drops  sizes 
and  allow  input  for  wind,  which  is  represented  by  a  3-dimesional  vector  applied  to  every 
drop.  The  wind  vector  is  implemented  with  a  Gaussian  spread  among  drops  enabling 
investigation  of  correlation  effects,  which  can  be  applied  within  or  between  pulses. 

The  results  from  running  the  simulator  show  that  using  the  model  developed  will 
generate  data  following  expected  theoretical  distributions.  The  simulator  treats  the  rain 
drops  as  point  targets  adding  the  electric  fields  to  generate  the  RCS.  The  results  verify 
the  assumption  that  the  average  return  from  rain  can  be  approximated  by  adding  the 
RCSs  of  all  drops.  The  simulator  can  be  adapted  to  explore  effects  of  changes  in  DSD. 
Since  the  simulator  distributes  the  drops  randomly  within  the  whole  volume,  including 
the  resolution  volume  and  its  margins,  the  average  distribution  for  sub-cells  does  not  fully 
follow  that  of  Marshal-Palmer.  The  effects  of  this  should  be  further  explored. 

The  current  version  of  the  simulator  does  not  take  into  account  attenuation, 
mutual  coupling,  and  multiple  scattering.  Neither  losses  nor  noise  are  introduced  and 
dual  polarization  is  left  for  future  work.  Although  simplified,  the  weather  radar  signal 
simulator  based  on  a  physical  model  of  a  spatial  region  containing  raindrops  produced 
from  this  research  work  will  allow  for  further  development  improving  realism. 
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B.  CAN  REALISTIC  MOTION  THAT  DEPENDS  UPON  THE  TYPE  OF 

WEATHER  SYSTEM  BE  IMPARTED  TO  THE  RAINDROPS? 

This  work  has  prioritized  the  electromagnetic  and  radar  aspects  and  the 
meteorology  portions  have  been  kept  to  a  minimum.  The  simulator  version  presented 
uses  a  vector  model  to  introduce  motion  to  the  raindrops.  This  can  be  adapted  to  impart 
realistic  motion.  To  do  that,  realistic  models  need  to  be  acquired  in  vector  form  allowing 
implementation  in  the  simulator.  The  current  version  of  the  simulator  measures  only  one 
resolution  cell  exploring  effects  of  de-correlation  using  different  schemes.  However,  the 
simulator  can  be  extended  to  cover  a  larger  volume  including  several  range  bins.  Using 
velocity  spread  as  an  input  parameter  makes  it  possible  to  simulate  spectrum  widening 
effects  without  having  to  create  a  detailed  physical  model  of  each  spreading  factor. 
There  is  no  doubt  that  using  a  vector  model  would  make  it  possible  to  implement  realistic 
motion  to  the  raindrops  although  this  has  not  been  done  in  this  version  of  the  simulator.  It 
was  beyond  the  scope  of  the  thesis  to  fully  implement  realistic  motion  models. 


C.  CAN  THE  OUTPUT  OF  THE  WEATHER  RADAR  SIMULATOR  BE 
USED  TO  STUDY  SPATIAL  AND  TEMPORAL  SAMPLING  SCHEMES 
FOR  CONSTANT  OR  STEPPED  FREQUENCY  SAMPLING  PULSES? 

Two  methods  have  been  implemented  focusing  on  temporal  sampling;  namely 

adapting  Pulse  Repetition  Frequency  so  the  drops  have  time  to  reshuffle  between  pulses, 

and  Frequency  shifting  between  sample  pairs.  A  simulation  was  run  with  the  following 

data  (for  details  see  appendix  C): 


• 

Wavelength,  A 

0.1  m 

• 

PRF 

300  Hz 

• 

Wind 

-5  m/s  in  x-direction  (=radial  velocity) 

• 

Velocity  spread 

2  m/s 

To 

determine  a  suitable  PRF  the  sample  time  threshold  provided  in  Equation 

(2.18)  was  used.  The  average  Z  of  the  simulator  data  lies  within  0.5  dB  of  the  integral 
form  which  is  also  true  for  the  value  retrieved  using  the  power  samples  (see  Figure  25). 
The  reflectivity  estimate  based  on  sample  power  lies  approximately  1  dB  below  the 
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estimated  average  reflectivity,  Z  and  0.5  dB  below  the  ground  truth,  which  can  be 
explained  by  the  approximation  of  the  radar  beam  where  only  drops  inside  the  half  power 
beamwidth  are  considered.  The  drop  position  within  the  radar  beam  will  also  affect  the 
gain  level,  which  makes  the  position  of  the  drops  vital.  A  constant  gain  approximation 
might  result  in  a  more  accurate  approximation.  The  chosen  PRF  led  to  a  level  of 
correlation  of  p~  0.543,  which  is  just  below  the  threshold  of  0.6  but  still  enough  to 
make  it  difficult  to  estimate  radial  velocity  with  reasonable  accuracy.  The  sample  pairs 
are  not  strongly  correlated  which  makes  an  estimate  possible  but  the  variance  will  cause 
errors  partly  due  to  the  lack  of  correlation,  and  partly  because  there  are  too  few  samples. 

A  new  simulation  was  run  this  time  using  a  rapid  scanning  scheme  with  stepped 
frequencies  as  described  in  a  previous  section.  The  new  simulation  was  run  using  the 
same  inputs  but  adding  another  PRF  of  5000  Hz  to  improve  the  Doppler  estimation 
capability. 

The  result  were  very  similar  to  the  first  run  using  only  time  to  acquire  de- 
correlation,  however,  comparing  the  total  time  to  produce  50  samples  differ  significantly, 
still  accomplishing  a  Z-value  close  to  the  integral  form  (Figure  26).  The  reason  for  this  is 
that  the  300  Hz  PRF  only  needs  to  be  used  to  separate  equal  frequencies  while  stepping 
between  frequencies  can  be  done  very  quickly.  Three  de-correlated  samples  can  now  be 
collected  in  almost  the  same  amount  of  time  as  previously  only  one  pulse  could  be.  The 
pulse  pair  scheme  also  enables  an  estimation  of  radial  velocity  seen  in  Figure  27.  The 
spread  was  larger  than  the  input  of  2  m/s  which  can  be  explained  by  the  fact  that  the 
spread  was  applied  both  between  pulses  and  within  pulses  between  different  drop  sizes. 
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Estimated  Z  per  pulse 
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Figure  25.  Temporal  Sampling,  PRF  variation.  PRF  300  Hz.  Velocity 
Spread  2  m/s  between  pulses  and  within  pulses  for  drops  of  different  sizes. 


The  output  from  the  latest  simulation  presented  an  average  radial  velocity  of 
-5.25  m/s  and  a  velocity  spread  of  4.6  m/s  (using  Equations  (2.33)  and  (2.34))  which  is 
close  to  the  expected  values. 
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0  10  20  30  40  50  60 

Sample  time  [ms] 


Figure  26.  Temporal  Sampling,  Stepped  Frequency.  PRFi  5000  FIz,  PRF2 
300  FIz.  Velocity  Spread  2  m/s  between  pulses  and  within  pulses  for  drops  of 

different  sizes. 
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Figure  27.  Velocity  histogram  with  Gaussian  Fit,  Temporal  Sampling 
using  stepped  frequency.  The  input  velocity  was  -5  m/s.  dv  =  4.6  is  estimated 

using  Equations  (2.33)  and  (2.34). 
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To  analyze  the  simulator  behavior  over  time  for  several  consecutive  runs,  a 
simulation  series  was  performed  using  50  runs  and  the  results  were  analyzed  for  different 
number  of  pulses.  In  [2]  the  intensity  probability  distribution  for  the  averaged  return  for 
different  numbers  of  pulses  was  plotted  as  in  Figure  28,  where  the  dependence  of  the 
number  of  pulses  can  be  seen. 


Figure  28.  Probability  distribution  of  Jk,  which  is  the  return  intensity. 
The  intensity  average  is  A 2  (From  [2]). 

Running  the  simulator  using  a  A  =  0.1  m,  and  a  PRF  of  200  Hz,  which  represents 
a  correlation  coefficient  of  about  0.04  (Equation  (2.15))  considering  a  velocity  spread  of 
4  m/s,  renders  results  seen  in  Figure  29.  The  figures  resemble  each  other  in  shape  but  for 
A  =  100  samples,  the  peak  for  the  simulator  distribution  is  not  as  high  as  the  result  in 
Figure  28.  As  the  number  of  pulses  goes  up  the  distribution  become  more  Gaussian, 
which  is  expected  for  any  distribution  if  the  number  of  samples  increases.  Also,  the 
spread  decreases  as  the  number  of  pulses  increases.  Having  completely  independent 

samples  will  change  the  average  velocity  spread  as  <7-  =  ~^= ,  where  k  is  the  number  of 
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pulses  or  samples.  The  results  are  somewhat  expected  since  the  distribution  of  the  power 
return  is  expected  to  be  exponential  [3].  The  standard  deviation  of  an  exponential 
distribution  is  expected  to  be  equal  to  the  mean  of  the  same,  which  the  results  verily. 

However,  if  samples  are  independent  the  spread  should  decrease  with  a  factor  of  1  /  \[k 

as  k  increases.  For  2  pulses  the  expected  value  equals  1/  v2  =  0.71 ,  which  matches  the 
simulator  result.  For  10  pulses  the  value  is  0.316,  while  the  simulator  shows  0.41. 
Finally  for  100  pulses  the  value  is  0.10.  The  simulator  shows  0.21,  which  is  more  than 
the  double  what  is  expected.  This  suggests  that  not  all  samples  can  be  considered 
independent.  Why  this  occurs  when  a  sufficient  spread  is  applied  cannot  be  explained  at 
this  point  but  is  of  relevance  for  future  investigations. 

Histogram 


1  pulse 

-  —2  pulses 

-  -  10  pulses 

-  -  -100  pulses 


Figure  29.  Histogram  showing  distribution  of  intensity  (power)  from  200 
consecutive  runs.  Number  of  pulses  averaged  is  1,  2, 10  and  100.  p  ~  0.04  so 

independent  samples  are  expected.  The  relative  standard  deviation  (<x  r  IP  )  in 
each  case  is,  1  pulse=0.99,  2  pulses=0.74, 10  pulses=0.41,  and  100  pulses=0.21. 
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D.  CAN  THE  WEATHER  RADAR  SIMULATOR  BE  USED  TO  STUDY  THE 
USE  OF  PULSE  COMPRESSION  AND  RANGE  AVERAGING  AS  A 
MEANS  OF  RAPIDLY  OBTAINING  INDEPENDENT  SAMPLES? 

A  simulation  was  performed  using  the  simulator  version  adapted  to  pulse 

compression  as  described  in  a  previous  section.  The  resolution  volume  differs  from  the 

previous  simulations  because  of  the  stretched  pulse  making  the  range  dimension  of  the 

resolution  cell  10  times  larger.  All  other  data  was  set  as  before,  i.e.  PRF  of  300  Hz,  wind 

speed  of  5  m/s  and  a  velocity  spread  of  2  m/s.  The  number  of  pulses  is  still  50  but  they 

are  used  differently  letting  every  pulse  cover  a  range  bin  equal  to  previous  runs.  The  new 

stretched  pulse  consists  of  10  sub-pulses  of  different  frequencies  creating  the  chirp.  Each 

new  pulse  is  summed  over  all  frequencies  stepped  over  the  whole  resolution  cell  and  then 

averaged.  Effectively  only  5  pulses  are  used,  but  since  each  pulse  consists  of  10 

frequency  steps,  the  average  value  is  based  on  50  samples. 


Estimated  Z  per  pulse 


Figure  30.  Pulse  Compression  simulation.  5  pulses,  Frequency  stepped 
from  3  GHz  to  3.0018  GHz  with  0.2  MHz  steps  for  each  pulse. 

The  estimate  of  Z  based  on  averaging  lies  1.0  dB  above  the  result  using  the 
integral  form,  which  indicates  good  accuracy.  The  smaller  deviation  between  samples 
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compared  with  previous  simulations  can  be  explained  by  the  fact  that  each  value  really 
represents  an  average  of  10  range  bins.  The  higher  level  can  be  explained  by  the  method 
used  to  average.  The  volume  of  the  different  resolution  bins  increases  moving  away  from 
the  radar  so  if  this  is  not  taken  into  account  when  averaging,  the  values  will  end  up 
deviated  from  the  expected  value.  Another  explanation  is  the  fact  that  the  spread  in 
frequency  is  not  enough  to  ensure  sufficient  de-correlation  to  get  independent  samples. 
Running  the  simulator  for  8  stretched  pulses  generated  a  result  very  close  to  the  integral 
form  (0.4  dB)  something  that  needs  further  investigation  for  verification. 

The  results  show  the  possibility  to  use  the  simulator  to  study  the  use  of  pulse 
compression  and  range  averaging  as  a  means  of  rapidly  obtaining  independent  samples. 


E.  CAN  THE  OUTPUT  OF  THE  WEATHER  RADAR  SIGNAL  SIMULATOR 

BE  USED  TO  PRODUCE  DISPLAYS  OF  PEDAGOGICAL  INTEREST? 

During  the  progress  of  building  the  simulator,  several  displays  have  been 
developed  and  used  to  assist  with  the  evaluation  of  each  module.  Once  the  simulator 
output  data  has  been  produced  it  is  all  numbers  making  it  possible  to  post  process  in 
many  different  ways,  displays  being  one.  To  be  able  to  use  the  simulator  in  teaching, 
visual  aids  presenting  the  signal  pulse-to-pulse  were  developed  as  well  a  phasor 
representation  displaying  the  phase  change  over  time  between  pulse  pairs.  Adding  a 
Graphical  User  Interface,  GUI,  would  probably  enhance  the  utility  the  simulator  as  a 
pedagogical  tool  for  investigating  the  effects  of  different  weather  phenomena  and 
variation  of  input  values  for  both  radar  and  weather  parameters.  Nevertheless,  the  output 
of  the  existing  simulator  can  be  used  to  produce  displays  of  pedagogical  interest. 


F.  CAN  THE  OUTPUT  OF  THE  WEATHER  RADAR  SIGNAL  SIMULATOR 
BE  USED  TO  STUDY  THE  PERFORMANCE  OF  ESTIMATORS  FOR 
THE  FIRST  THREE  DOPPLER  MOMENTS? 

The  first  versions  of  the  simulator  only  allowed  for  changes  in  a  single  PRF 
making  it  possible  to  vary  time  between  pulses  effecting  de-correlations.  The  simple 
form  of  return  signal  produced  samples  equally  spaced  in  time  making  it  simple  to 
perform  post  processing  in  the  form  of  auto  correlation  and  FFT.  The  spectrum  was 
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plotted  and  it  was  possible  to  derive  an  estimated  radial  velocity  of  the  drops  from  the 
average  Doppler  frequency.  When  the  complexity  of  the  return  signal  increased  and  the 
simulator  allowed  for  multiple  PRFs  and  delays,  the  signal  needed  manipulation  to  enable 
FFT  or  autocorrelation.  Initial  trials  padding  the  signal  with  zeros  making  the  samples 
equally  spaced  in  time  were  performed  but  adding  time  to  a  signal  demands  more  signal 
processing  to  remove  effects  of  folding  and  aliases.  Digital  signal  processing  is  beyond 
the  original  scope  of  the  thesis  which  is  why  this  interesting  trail  was  abandoned.  There 
are  no  test  runs  verifying  that  the  simulator  can  be  used  to  study  the  performance  of 
estimators  for  the  first  three  Doppler  moments  but  that  does  not  imply  that  the  simulator 
cannot  be  used  for  this  application.  Having  a  complex  return  signal  allows  for  post 
processing  of  many  sorts,  as  mentioned  earlier. 


G.  CAN  THE  OUTPUT  OF  THE  WEATHER  RADAR  SIGNAL  SIMULATOR 
BE  USED  TO  STUDY  THE  UTILITY  OF  ESTIMATORS  FOR  HIGHER 
ORDER  DOPPLER  MOMENTS? 

Since  the  estimation  of  the  lower  order  of  estimators  was  not  performed,  the 
higher  order  estimators  have  not  been  investigated. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE 

DEVELOPMENT  AND  WORK 


A.  CONCLUSIONS 

The  objective  of  the  research  was  to  build  a  weather  radar  signal  simulator 
delivering  I  and  Q  channel  outputs  from  a  physical  representation  of  rain  drops.  A 
previously  developed  simulator  was  initially  evaluated  and  verified  serving  as  the  basic 
building  block  when  developing  the  final  version.  The  Radiation  Integrals  were  used  to 
derive  the  electromagnetic  scattering  from  rain  drops  approximated  as  dielectric  spheres. 

Without  accounting  for  multiple  scattering,  attenuation,  or  coupling  effects,  the 
scattered  electric  field  from  each  drop  is  summed  for  every  pulse  providing  the  total  RCS 
of  the  rain  in  the  radar  resolution  cell  as  well  as  a  complex  voltage  return  representation. 
Using  the  Marshall-Palmer  drop  size  distribution,  an  initial  position  for  every  drop  is 
determined  within  an  extended  volume  enabling  the  drops  to  fall  through  the  radar  beam 
for  all  samples.  Between  pulses  the  drops  are  moved  according  to  their  respective 
tenninal  velocities  based  on  the  Atlas-Ulbrich  approximation  and  wind.  To  allow  for 
width  in  the  velocity  spectrum,  a  Gaussian  velocity  spread  is  applied  to  the  drops 
between  and  within  pulses.  The  simulator  power  return  is  weighted  using  the  Radar 
Range  Equation  for  point  targets,  considering  all  drops  as  point  targets,  and  evaluated  as 
a  volume  target  using  the  Weather  Radar  Range  Equation. 

Several  different  functionalities  have  been  implemented  allowing  for  stepped 
frequencies,  multiple  PRFs,  pulse  compression  using  a  chirp,  and  variation  of  radar  and 
weather  input  parameters.  Post  processing  capabilities  include  autocorrelation  and  FFT 
(only  for  single  PRF  version);  computation  of  weather  parameter  estimates  such  as 
reflectivity  factor,  Z;  average  doppler,  radial  velocity,  and  velocity  spread;  pedagogical 
plots  including  a  Phasor  plot  of  phase  change  over  time,  a  velocity  histogram,  power  and 
reflectivity  for  each  pulse  over  time. 

This  research  has  shown  that  building  a  weather  radar  signal  simulator  based  on  a 
physical  model  of  a  spatial  region  containing  raindrops  is  possible.  The  benefit  of  using  a 
physical  model  is  improved  realism  and  the  ability  to  study  effects  of  scientific  interest. 
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Also,  the  model  permits  isolation  of  effects  making  it  possible  to  explore  the  results  when 
varying  some  parameters  while  keeping  others  fixed.  Using  a  vector  model  to  input  wind 
and  other  effects  on  the  rain  drops  makes  it  possible  to  impart  realistic  motion  to  the 
drops  although  this  has  not  been  done  fully  in  this  version  of  the  simulator.  This  was 
beyond  the  scope  of  this  thesis  but  is  a  recommended  area  for  further  study.  Evaluating 
the  result  from  averaging  radial  velocity  rendered  an  interesting  finding  namely  the 
difference  between  using  FFT,  autocorrelation,  and  pulse  pair  processing  to  estimate 
radial  velocity.  Using  the  FFT  or  autocorrelation  skewed  the  result  towards  the  bigger 
drops  because  the  amplitude  infonnation  is  kept  in  the  process.  On  the  other  hand  with 
the  pulse  pair  algorithm,  only  the  phase  information  is  used,  resulting  in  a  “true”  average 
velocity  for  all  drops  regardless  of  size.  Intuitively  smaller  drops  would  be  more 
susceptive  to  wind  than  bigger  drops  which  would  emphasize  the  need  for  pulse  pair 
algorithms  over  the  FFT  or  autocorrelation  for  better  wind  velocity  estimates.  No 
theoretical  investigation  or  look  for  research  in  this  area  has  been  perfonned  but  is 
suggested  for  future  studies.  Different  implementations  of  velocity  spread  have  been 
tested  leading  to  variations  in  results.  The  final  version  of  the  simulator  allows  for  spread 
both  between  pulses,  which  represents  how  the  average  velocity  is  spread  between 
samples;  and  within  pulses  with  a  size  dependent  spread,  which  represents  the  spread 
within  the  resolution  cell. 

The  title  of  this  research  indicates  where  the  priority  lies,  namely  developing  the 
simulator  to  examine  sampling  rates  and  scanning  schemes.  The  results  show  that  it  is 
possible  to  use  the  weather  radar  simulator  to  study  spatial  and  temporal  sampling 
schemes  for  constant  or  stepped  frequency  sampling  pulses.  The  research  aimed  at 
building  the  radar  and  not  to  fully  explore  the  effects  of  different  rate  and  schemes,  but 
the  results  from  testing  and  verification  indicate  that  using  methods  to  more  rapidly 
obtain  independent  samples  decreases  the  time  to  obtain  accurate  estimates  of  weather 
signal  parameters.  By  introducing  three  frequency  steps,  the  time  to  acquire  accurate 
data  was  reduced  by  a  third,  still  being  able  to  acquire  Doppler  information  as  well  as 
reflectivity.  However,  more  simulations  need  to  be  performed  comparing  results  when 
parameters  are  varied  to  fully  explore  the  benefits  of  frequency  stepping  and  its 
limitations. 
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By  implementing  a  pulse  compression  capability,  using  frequency  stepping  and 
range  averaging,  results  were  provided  showing  the  ability  to  rapidly  obtain  independent 
samples  using  pulse  compression.  Increasing  the  radar  resolution  cell  slows  down  the 
simulator,  which  can  be  compensated  by  increasing  the  Norm- factor  (see  previous 
chapters  for  details).  Again  further  testing  is  suggested  exploring  the  effects  of  varying 
parameters. 

Post  processing  the  signal  samples  provides  wide-ranging  opportunities  for 
different  displays  of  pedagogical  interest.  Most  plots  have  been  developed  to  compare 
and  verify  results  but  some  plots  have  been  developed  only  to  display  pedagogically 
interesting  phenomena.  Being  able  to  follow  phase  change  over  time  is  an  example  of  the 
latter.  Implementing  a  GUI  will  probably  be  the  greatest  improvement  in  this  area. 

Having  an  accurate  signal  allows  for  more  than  just  plots  of  results,  it  opens  up 
possibilities  for  post  processing  schemes  beyond  the  scope  of  this  work.  Previously 
indicated  processing  includes  filtering,  FFT  and  autocorrelation.  In  the  introduction, 
cross  scientific  benefits  were  mentioned  signifying  the  use  of  knowledge  from  one  area  in 
another.  Having  a  weather  signal  consisting  of  rain  echoes  can  be  used  to  show  or 
investigate  filtering  to  remove  weather  clutter. 

The  initial  intentions  included  ideas  of  improving  the  simulator  speed  by 
optimizing  algorithms  and  scripts.  This  has  not  been  accomplished.  However,  the 
simulator  runs  a  simulation  locally  on  a  PC  with  an  Intel  Pentium  4  processor  550 
(3.4GHz,  1M,  800MHz  FSB,  NTFS  Memory  4.0GB  DDR2  Non-ECC  SDRAM,  400Hz 
(4DIMM)),  in  just  a  few  minutes,  depending  on  complexity  and  input  parameters. 


B.  RECOMMENDATIONS  FOR  FUTURE  RESEARCH 

The  simulator  is  ready  to  do  its  job,  setting  more  than  one  path  for  future  work. 
The  focus  in  this  research  has  been  to  develop  a  tool  for  research,  which  leaves  the 
weather  research  parts  for  future  work.  Two  main  areas  of  development  and  research  are 
proposed: 

•  Examine  sampling  rates  and  scanning  schemes  for  fast  volumetric  scanning 
and  update 
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Simulator  development 


1.  Examine  Sampling  Rates  and  Scanning  Schemes  for  Fast  Volumetric 

Scanning  and  Update 

Weather  Radar  Signal  Simulator  should  be  run  to  examine  weather  signal 
parameter  estimation  sensitivity  when  varying  radar  parameters  and  precipitation 
estimator  models  where  changing  the  drop  size  distribution  would  be  of  interest. 
Comparison  should  be  made  to  real  life  measures  as  well  as  theoretical  models.  To 
complement  this  work  weather  models  should  be  explored  to  find  the  physics  behind 
different  weather  phenomena  enabling  implementation  in  the  simulator.  Alterations  in 
the  simulator  might  be  necessary  to  facilitate  this  implantation. 

2.  Simulator  Development 

Even  though  the  final  version  of  the  simulator  has  several  interesting 
functionalities  many  interesting  effects  and  functionalities  can  be  added. 

a.  Adding  a  GUI 

The  main  priority  developing  the  simulator  should  be  to  create  a  GUI  making  the 
tool  more  user-friendly. 

b.  Modeling  Non-Spherical  Drops 

The  current  simulator  approximates  each  drop,  regardless  of  size,  as  a 
dielectric  sphere.  Research  show  that  drops  flatten  [9]  while  falling  so  an  oblate  spheroid 
approximation  would  be  more  appropriate.  Work  done  by  Bringi  and  Chandrasekar  [14] 
provide  good  mathematical  tools  enabling  implementation  in  the  simulator. 


c.  Adding  Dual  Polarization 

Implementing  a  spheroid  approximation  makes  dual  polarization 
measurements  possible  and  new  estimates  can  be  computed  such  as  differential 
reflectivity.  The  complexity  of  this  can  be  overwhelming  considering  effects  like  cross 
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talk  and  multiple  scattering.  This  is  a  problem  for  scientists  or  students  who  are 
searching  for  a  real  challenge. 

d.  Post  Processing  for  Higher  Order  Moments 

The  output  signal  can  be  post  processed  to  produce  the  power  spectrum, 
which  was  done  in  the  early  versions  that  had  uniform  time  spacing  between  samples. 
Applying  digital  signal  processing  to  the  output  signal,  retrieving  the  power  spectrum 
would  make  it  possible  to  explore  higher  order  moments  of  the  spectrum  and  to  evaluate 
their  utility  as  estimators. 

e.  Linking  Simulator  Input  to  Storm  Physics 

As  mentioned,  stonn  models  should  be  explored  to  detennine  the  physics 
associated  with  different  weather  phenomena.  The  physics  can  be  used  to  detennine  the 
wind  velocity  and  velocity  spread  that  should  be  used  as  simulator  input  parameters. 
This  will  add  to  the  realism  of  the  simulator  and  will  make  it  possible  to  explore  effects 
on  the  radar  signal  caused  by  different  weather  phenomenon. 

f  Adding  Noise  to  the  Weather  Signal 

Real  weather  signals  will  always  be  contaminated  by  some  noise.  Noise 
from  all  sources  should  therefore  be  added  as  a  set  of  inputs.  The  benefit  of  this  would 
be  the  ability  to  study  the  degradation  caused  by  noise. 

g.  Adding  Mutual  Coupling 

Mutual  coupling  was  examined  but  needs  more  attention  and  different 
techniques  to  enable  implementation  into  the  simulator.  For  large  rainfall  rates,  raindrops 
may  be  spaced  at  distances  that  result  in  mutual  coupling.  The  total  RCS  in  this  case  will 
be  less  than  that  computed  assuming  no  coupling.  This  effect  needs  further  investigation 
to  determine  what  stonn  types  and  rain  rates  result  in  coupling  that  renders  the 
independents  scattering  assumption  invalid. 
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APPENDIX  A.  DERIVATION  OF  RCS  FOR  SINGLE  DROP 


A  water  drop  can  be  approximated  as  a  small  dielectric  sphere.  Let  the  sphere 
radius  be  a,  and  the  permittivity  e ,  and  penneability  z/0 ,  also  let  the  sphere  be  centered  at 
the  origin.  An  incident  plane  wave  is  present 

E,  =  zE0e~jkx 

The  electric  field  at  a  point  in  space  can  be  expressed  as 

E  =  Ei+  Es 

Using  the  volume  equivalent  principal,  the  current  J  can  be  written  as 

J  =  jco(e-£0)E 

Letting  a » X ,  E  is  then  approximately  E  =  zEl  (constant  inside).  Using 

Es  =  —  E]  can  be  derived  as 
3  CO£0 


zE\  -  E  =  Ei  +  Es  =  Ei  + 


JJ 

3(D£n 


zEi  =  Ei  +  j 


jCQ(£,.£0-£0)E 


0  °'--Ei- 


(£,.  —  X)zE\ 


3co£„ 


zE  i  = 


3  Ei 

(£,.+2) 


Now  solving  for  J 


zEi 


3  Ei 

(£,.  +  2) 


jO)(£~£0 ) 

j_j3CQ£0(£r-l)Ei 

(£,.+2) 
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Finally  putting  the  current  in  the  radiation  integral 


F.(r,6A)  =  Je^dv,  g  =  0 


E(r,e,t)  =  =Ae-* 


4;rr 


4;rr 

_7,T  j3(O£0(£r  -\)Ei  | ^ 


a  a  a 


(£r  +  2) 


0  0  0 
3 


2(r,e,*)  =  liHe-J*  j^^-DE.^a 


4nr 


(£r  +  2) 


E(r,0, (/))  =  -  [Ke-J^  2^c£o(£r  y)E‘ 


r  V  £„ 


A(sr  +  2) 


E(r,0,(/))  =  —e-Jkr^—^Eia3 
r  (er+2) 


Using  this  in  the  formula  for  RCS  renders,  (letting  a  —  D/2) 


<Jh  =  4  nr 


Ei 


=  4  nr2 


k2  c_jkr  (er  -1) 

(U+2) 


Eta 


4nk  a  (£r-l)~ 

(£,  +  2)2 


^■5(f,.-l)2 


/l4  (£•,.  +  2)2 


-D 
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APPENDIX  B.  MICROWAVE  STUDIO  SIMULATION  SETUP 


First  the  target(s)  was  created.  All  targets  were  spheres  of  4  mm  in  diameter  with 
the  specifications  of  a  water  drop.  The  incident  ray  was  a  plane  wave  in  the  z-direction 
with  polarization  in  the  x-direction  or  in  the  y-direction.  When  adding  the  second  drop, 
the  drop  was  placed  side  by  side  with  the  first  drop,  i.e.  in  the  x-direction.  The  plots  are 
farfield  plots  in  spherical  coordinates.  The  mesh  resolution  was  set  to  20  lines  per 
wavelength. 


Figure  1.  Target.  Water  drop. 
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Figure  2.  Incident  Plane  for  the  one-drop  simulation 


Figure  3.  Type  of  plot  and  definition  of  angles. 
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Hc^lipluiie  ut  V«  9  (  I ndux=  8  ) 


Figure  6.  Resolution  of  the  Meshed  grid  for  two  drops. 


Figure  7.  Incident  Plane  for  the  two-drop  simulation 
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Figure  9.  Plot  of  results  adding  extra  volume  when  calculating  RCS  for  one  drop  in  the  0  - 

plane  for  (j)  =  0 . 
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APPENDIX  C. 


Radar  Data 

P_t=10e3; 

G_dB=20; 

G=10A(G_dB/10); 

angle_3db=0.20; 

pulsewidth=0.20; 

elevation=3; 

margine=5; 

azimuth=0; 

number_pulses=5  0 ; 

PRT=l/300; 

distance=0.5; 

wavelength=0.1; 

k=2*pi/wavelength; 

Rain  parameters 

diam_low=0.05; 

diam_lim=l; 

diam_mid=5; 

diam_high=T5; 

resolutionl=0.08; 

resolution2=0.3; 

resolution3=0.8; 

sort_approx=l; 

kw=0.93; 

rain_rate=50; 

rain_anglel=0; 

rain_angle2=0; 

wind_vector=[-5,0,0]; 

spread=2; 

Norm=100; 


SIMULATION  SETUP  FOR  THE  WEATHER 
RADAR  SIMULATOR 


%  Transmitted  power 
%  Maximum  gain  in  dB 
%  Maximum  gain 

%  3  dB  Beam  angle  [degrees],  (half  half  power  beamwidth) 
%  Pulsewidth  [microsecond] 

%  Elevation  angle  [degrees] 

%  margin  for  resolution  cell[m] 

%  Azimuth  angle  [rad] 

%  Number  of  pulses 
%  Pulse  Repetition  Time  [seconds] 

%  Range  from  the  radar  to  the  resolution  cell  [Km] 

%  Wavelength  [m] 


%  Lowest  diameter  of  Precipitation  [mm] 

%  Divider  between  low  and  mid  [mm] 

%  Divider  between  mid  and  high  [mm] 

%  Highest  diameter  of  Precipitation  [mm] 

%  Lower  diameter  resolution  [mm] 

%  Higher  diameter  resolution  [mm] 

%  Highest  diameter  resolution  [mm] 

%  RCS  approximation  model;  [l]=Rayleigh 
%  abs(K)A2  approx  for  water 
%  Rain  Rate  [mm/h] 

%  degrees  Rain  Angle  (Phi) 

%  degrees  Rain  Angle  (Theta) 

%  Wind  speed  vector  (x,y,z)  [m/s] 

%  Wind  speed  spread  [m/s],  (all  directions) 
%  Normalization  coefficient 
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APPENDIX  D.  CODE  DEVELOPED  IN  MATLAB. 


There  are  currently  four  versions  of  the  simulator: 


1)  The  basic  version,  WR  basic,  that  will  run  with  one  PRF.  The  plots  include 

•  Doppler  per  pulse  pair  with  average  Doppler, 

•  Radial  velocity  histogram  with  Gaussian  fit 

•  Return  power  with  average  Return  power 

•  Reflectivity  estimations  with  Ground  Truth  for  comparison 

•  Voltage  Correlation  function  with  estimation  of  p[r  =  TS) 

•  Power  spectrum  of  return  signal 

2)  Dual  PRF  version,  WR_dual_prf,  that  will  run  with  two  PRFs.  The  plots 

include 

•  Doppler  per  pulse  pair  with  average  Doppler, 

•  Radial  velocity  histogram  with  Gaussian  fit 

•  Return  power  with  average  Return  power 

•  Reflectivity  estimations  with  Ground  Truth  for  comparison 

•  Phasor  plots  for  the  30  first  pulse  pairs. 

3)  Frequency  step  version,  WR_fq_step,  that  will  step  through  several  frequencies 
and  transmit  a  pulse  pair  on  every  frequency.  One  PRF  controls  the  time  between  pulses 
in  the  pulse  pair  and  one  PRF  controls  the  time  between  equal  frequencies.  The  plots 
include 

•  Doppler  per  pulse  pair  with  average  Doppler, 
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•  Radial  velocity  histogram  with  Gaussian  fit 

•  Return  power  with  average  Return  power 

•  Reflectivity  estimations  with  Ground  Truth  for  comparison 

•  Phasor  plots  for  the  30  first  pulse  pairs. 


4)  Pulse  compression  version,  WR_pulse_com,  that  divide  the  resolution  volume  into 
range  bins  and  averages  over  bins  and  pulses.  The  chirped  pulse  is  approximated  by 
discrete  stepped  frequencies.  The  plots  include 

•  Return  power  with  average  Return  power 

•  Reflectivity  estimations  with  Ground  Truth  for  comparison 

1)  WR_basic 


%  Weather  Radar  Signal  Simulator,  developed  version 
%  Ulf  Schroder 
%  July,  2005 

% - 


%  Main  program  for  The  NPS/Marielle  Gosset  Weather  Radar  Signal  Simulator. 


% - 

%  DESCRIPTION 
%  Part  I 

%  From  Radar  Data: 

%  -  Generate  Beam  Resolution  cell  approximation 

%  -  Generate  Cube  adding  a  margine  to  the  Beam  Resolution  cell  to  allow  for 
%  the  rain  drops  to  fall  trough  the  through  beam  over  time. 

% 

%  Part  II 

%  From  Rain  Parameters: 

%  Generating  the  fundamental  T  Matrix  concisting  of 
%  -  Diameter  Vector  and  Number  of  drops  in  each  diameter  interval 
%  using  the  Marshall-Palmer  exponential  approximation. 

%  -  the  backscatter  RCSs  for  each  Diameter  size 
%  -  the  total  number  of  drops  for  each  size  per  mA3 

%  Part  III 

%  Randomly  place  all  drops  in  cube 
%  Part  IV 

%  For  every  pulse  the  coherent  Electric  field  return  is  calculated  Between 
%  pulses  all  drops  are  moved. 


%  Part  V 
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%  Calculate  and  plot  results 

% - 

clear  all 
close  all 
clc 

%  Input  Data 

c=3e8;  %  Speed  of  light  [m/s] 

%  Radar  Data 

P_t=10e3;  %  Transmitted  power 

G_dB=20;  %  Maximum  gain  in  dB 

G=10A(G_dB/10);  %  Maximum  gain 

angle_3db=0.20;  %  3  dB  Beam  angle  [degrees],  (half  half  power  beamwidth) 

pulsewidth=0.20;  %  Pulsewidth  [microsecond] 

elevation=3;  %  Elevation  angle  [degrees] 

margine=5;  %  margin  for  resolution  cell[m] 

azimuth=0;  %  Azimuth  angle  [rad] 

number_pulses=40;  %  Number  of  pulses 

PRT=l/200;  %  Pulse  Repetition  Time  [seconds] 

distance=0.5;  %  Range  from  the  radar  to  the  resolution  cell  [Km] 

wavelength=0.1;  %  Wavelength  [m] 

k=2*pi/wavelength; 

%  Rain  parameters 

diam_low=0.05;  %  Lowest  diameter  of  Precipitation  [mm] 
diam_lim=l;  %  Divider  between  low  and  mid  [mm] 

diam_mid=5;  %  Divider  between  mid  and  high  [mm] 

diam_high=15;  %  Highest  diameter  of  Precipitation  [mm] 

resolutionl=0.08;  %  Lower  diameter  resolution  [mm] 

resolution2=0.3;  %  Higher  diameter  resolution  [mm] 
resolution3=0.8;  %  Highest  diameter  resolution  [mm] 
sort_approx=l;  %  RCS  approximation  model;  [l]=Rayleigh 
kw=0.93;  %  abs(K)A2  approx  for  water 

rain_rate=50;  %  Rain  Rate  [mm/h] 

rain_anglel=0;  %  degrees  Rain  Angle  (Phi) 

rain_angle2=0;  %  degrees  Rain  Angle  (Theta) 
wind_vector=[-10;0;0];  %  Wind  speed  vector  (x,y,z)  [m/s] 
spread=l;  %  Wind  velocity  spread  [m/s],  (all  directions) 

spread_2=4;  %  Drop  velocity  spread  [m/s], 

Normal 00;  %  Normalization  coefficient 

%  Directional  cosines 

u=sin(pi/2-elevation*pi/180)*cos(azimuth*pi/180); 

v=sin(pi/2-elevation*pi/180)*sin(azimuth*pi/180); 

w=cos(pi/2-elevation*pi/180); 

% - 

%  Part  I,  Generate  the  resolution  volume 

% - 

%  Approximation  of  the  resolution  cell  of  the  pulse  circular  box, 

%  for  any  azimuth  angle 

% - 

%  INPUT 
%  -  distance 
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%  -  angle_3dB 
%  -  pulsewidth 
%  -  elevation 
%  -  margine 
%  -  azimuth 

% - 

%  OUTPUT 

%  -  volumeresolution:  computed  volume  of  the  resolution  cell 
%  -  volume_box:  volume  of  the  box  where  the  drops  are  put 
%  -  box:  the  8  rectangular  coordinates  of  the  box 
%  -  coord_vol_res:  the  coordinates  of  the  resolution  cell  where  the  drops 
%  are  taken. 

% - 


% - 

%  Creating  the  circular  beam  resolution  cell 
range_res=(c*pulsewidth*le-6)/2; 

%  Range  resolution  (ct/2)  [m] 

%  small  angle  approximation  allows  for  following 
angle_rad=(angle_3  db  *pi)/ 1 8  0 ; 

%  3  dB  beamwidth  [rad] 
radius  l=(tan(angle_rad)*(distance*  1000)); 

%  Radius  of  beam  at  R=distance 
radius2=(tan(angle_rad)*((distance*  1000)+range_res));  ... 

%  Radius  of  beam  at  R=distance  +  range  resolution 

surfacel=pi*radiuslA2; 

surface2=pi*radius2A2; 

%  volume  of  cone  =  base  surface  times  height  for  all  three 
volume_conel=(surfacel*(distance*1000))/3; 
volume_cone2=(surface2*((distance*1000)+range_res))/3; 
volume_resolution=(volume_cone2-volume_conel);  %  Beam  Cone  Volume 

%  computation  of  the  resolution  cell  coordinates,  spherical  coordinates 
%  elevation  in  radians 

elevation_rad=(elevation*pi)/180;  %  Elevation  [rad] 

azimuth_rad=(azimuth*pi)/180;  %  Azimuth  [rad] 

%  Matrix  used  when  comparing  wheather  a  drop  is  inside  or  outside 
%  beamwidth 

coord_vol_res=[distance*  1000, (distance*  1000)+range_res;elevation_rad-angle_rad,... 

elevation_rad+angle_rad;azimuth_rad-angle_rad,azimuth_rad+angle_rad]; 

% - 


% - 

%  Creating  the  box 

%  computation  of  the  box  and  its  coordinates 

coin  1 =[azimuth_rad+angle_rad,elevation_rad+angle_rad,  distance*  1 000] ; 

coin2=[azimuth_rad-angle_rad,elevation_rad+angle_rad,distance*1000]; 

coin3=[azimuth_rad+angle_rad,elevation_rad-angle_rad,distance*1000]; 

coin4=[azimuth_rad-angle_rad,elevation_rad-angle_rad,distance*1000]; 

coin5=[azimuth_rad+angle_rad,elevation_rad+angle_rad,distance*1000+range_res]; 

coin6=[azimuth_rad-angle_rad,elevation_rad+angle_rad,  distance*  1000+rangeres]; 

coin7=[azimuth_rad+angle_rad,elevation_rad-angle_rad,  distance*  1000+rangeres]; 
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coin8=[azimuth_rad-angle_rad,elevation_rad-angle_rad,distance*1000+range_res]; 


[xl,yl,zl]=sph2cart(coinl(l),coinl(2),coinl(3)); 

[x2,y2,z2]=sph2cart(coin2(l),coin2(2),coin2(3)); 

[x3,y3,z3]=sph2cart(coin3(l),coin3(2),coin3(3)); 

[x4,y4,z4]=sph2cart(coin4(l),coin4(2),coin4(3)); 

[x5,y5,z5]=sph2cart(coin5(l),coin5(2),coin5(3)); 

[x6,y6,z6]=sph2cart(coin6(l),coin6(2),coin6(3)); 

[x7,y7,z7]=sph2cart(coin7(l),coin7(2),coin7(3)); 

[x8,y8,z8]=sph2cart(coin8(l),coin8(2),coin8(3)); 

xmax=max([xl,x2,x3,x4,x5,x6,x7,x8]); 
xmin=min([xl,x2,x3,x4,x5,x6,x7,x8]); 
ymax=max([y  1  ,y2,y3,y4,y5,y6,y7  ,y8]); 
ymin=min([y  1  ,y2,y3,y4,y5,y6,y7  ,y8]); 
zmax=max([zl,z2,z3,z4,z5,z6,z7,z8]); 
zmin=min([zl,z2,z3,z4,z5,z6,z7,z8]); 


%  the  box  with  its  margin 

coin_l=[xmin-margine,ymin-margine,zmin-margine]; 

coin_2=[xmin-margine,ymax+margine,zmin-margine]; 

coin_3=[xmax+margine,ymax+margine,zmin-margine]; 

coin_4=[xmax+margine,ymin-margine,zmin-margine]; 

coin_5=[xmin-margine,ymin-margine,zmax+margine]; 

coin_6=[xmin-margine,ymax+margine,zmax+margine]; 

coin_7=[xmax+margine,ymax+margine,zmax+margine]; 

coin_8=[xmax+margine,ymin-margine,zmax+margine]; 


box=[coin_r,coin_2',coin_3',coin_4',coin_5',coin_6',coin_7',coin_8']; 

%  calculate  the  volume  of  the  box 

volume_box=abs((xmax+margine-(xmin-margine))*(ymax+8-(ymin-margine))*... 

(zmax+margine-(zmin-margine))); 


cube=[coin_r,coin_2',coin_3',coin_4',coin_r,coin_5',coin_6',coin_2',... 

coin_6',coin_7',coin_3',coin_7',coin_8',coin_4',coin_8',coin_5']; 

% - 


% - 

%  Part  II,  Generate  the  fundamental  T  Matrix 

% - 

%  This  part  produces  a  matrix  of  all  drops  that  contribute  to  the  total 
%  reflectivity  of  the  resolution  cell 

% - 

%  First  all  sizes  are  created,  with  two  vectors  of  different  size  resolution. 

%  The  size  vectors  represents  all  diameters.  Next  the  backscatter  RCS  for  each 
%  size  is  computed. 

%  Using  the  Marshall-Palmer  exponential  approximation,  the  number  of  drops  of 
%  each  diameter  is  computed. 

% - 

%  INPUT 
%  -  diam_low 
%  -  diam_lim 
%  -  diam_high 
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%  -  resolution  1 
%  -  resolution2 
%  -  sort_approx 
%  -  wavelength 
%  -  Kw 
%  -  rain_rate 

% - 

%  OUTPUT 
%  -  number_drops 
%  Matrix  T : 

%  -  first  row  all  sizes 
%  -  second  row  all  the  backscatter  RCSs 
%  -  third  row  the  total  number  of  drops  for  each  size  per  mA3 
%  the  columns  starting  from  the  smallest  size 

% - 

%  Create  dropsize  vector 

%  Precipitaion  Diameter  vector  (lowest<D<divider) 
vector_size_l=diam_low:resolutionl:diam_lim; 

%  Precipitaion  Diameter  vector  (divider<D<mid) 
vector_size_2=diam_lim:resolution2:diam_mid; 

%  Precipitaion  Diameter  vector  (mid<D<highest) 
vector_size_3=diam_mid:resolution3:diam_high; 

%  Precipitaion  Diameter  Vector 

vector_size=[vector_size_l,vector_size_2,vector_size_3]; 


% - 

%  Create  RCS  vector  for  all  Rain  drops  of  Precipitation  Vector. 
rcs=((kwA2)*(piA5)/(wavelengthA4))*((vector_size*le-3).A6); 
% - 


% - 

%  The  Marshall  Palmer  approximation  is  used  to  create  Precipitaion  Diameter 
%  Distribution  Vector  (Drop  Size  Distribution  N(D)  Vector) 

%  The  number  of  drops  are  scaled  down  by  a  factor  (Norm)  to  improve  the 
%  speed  of  the  simulator 

% - 

lambda=4.1*(rain_rateA(-0.21));  %[mm-l] 

number=8000*exp(-lambda*vector_size)/Norm; 

vector_resl=ones(l,length(vector_size_l))*resolutionl; 
vector_res2=ones(l,length(vector_size_2))*resolution2; 
vcctor_rcs3~oncs(  1 .  length)  vcctor_sizc_3))*reso  I  ution3; 

%  delta(Diameter)  dD 

vector_res=[vector_resl,vector_res2,vector_res3]; 

%  the  approximate  "true"  number  of  drops  N(D)dD  vector 
number=number.*vector_res; 

% - 


% - 

%  create  matrix  T 
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T=[vector_size;rcs;number]; 


% - 

%  Change  matrix  T  to  matrix  'matrixTRound' 

% - 

%  This  part  computes  the  total  number  of  drops  by  utilizing  the  resolution 
%  cell  and  the  matrix_T_Round 

%  The  only  difference  between  the  matrix_T_Round  and  T,  is  the  fact  that 
%  matrix_T_Round  has  the  total  number  of  drops  on  the  third  row  and  they 
%  are  rounded  of  to  nearest  integer  value. 

% - 


matrix_T_Round=T ; 

matrix_T_Round(3,:)=round(matrix_T_Round(3,:)*volume_box); 

%  Keep  track  of  errors  due  to  round  off 
diff=T(3 , :  )*volume_box-matrix_T_Round(3, : ); 

error_rcs=sum(diff.*matrix_T_Round(2,:))*Norm/volume_box*volume_resolution; 

error_Z=error_rcs/volume_resolution*wavelengthA4/piA5/abs(kw)A2/le-18; 


% - 

%  Part  III,  Initial  positioning  of  the  drops  within  the  volume 

% - 

%  This  part  generates  the  initial  position  of  the  drops  in  the  volume  of 
%  the  box. 

%  The  positions  are  recorded  in  the  file  "sizexx.maf  taking  into  account 
%  the  azimuth  angle 

% - 

%  INPUT 
%  -  box 

%  -  matrix_T_Round 
%  -  elevation 
%  -  azimuth 

% - 

%  OUTPUT 

%  All  files  including  the  sizes.  Each  contains  the  positions  of  all  drops 
%  for  a  given  diameter.  The  data  recoreded  as  follows: 

%  -  radial  distance  r 
%  -  elevation  phi 
%  -  azimuth  theta 

%  -  angle  with  respect  to  beam  center.  This  parameter  will  be 
%  used  to  identify  drops  which  are  in  the  beam. 

% - 


%  Randomly  put  the  calculated  number  of  drops  of  each  diameter  in  the  box 
number_size=length(matrix_T_Round(l,:)); 

for  i=l:number_size 

%  select  number  of  drops  of  certain  Diameter 
drop  s=matrix_T_Round(  3 ,  i ) ; 

%  create  as  many  positions  on  the  x  axis  as  #  drops  per  size 
elements_x=rand(  1  ,drops); 

%  put  them  in  the  range  of  interest 
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%  Random  placement  x  for  every  drop 

position_x=box(  1 , 1  )+(box(  1 ,3)-box(  1 , 1  ))*elements_x; 

%  create  as  many  positionson  the  y  axis  as  #  drops  per  size 
elements_y=rand(  1  ,drops); 

%  put  them  in  the  range  of  interest 
%  Random  placement  y  for  every  drop 

position_y=box(2, 1  )+(box(2,2)-box(2, 1  ))*elements_y; 

%  create  as  many  positions  on  the  z  axis  as  #  drops  per  size 
elements_z=rand(  1 ,  drops); 

%  put  them  in  the  range  of  interest 
%  Random  placement  z  for  every  drop 

position_z=box(3,l)+(box(3,5)-box(3,l))*elements_z; 

%  matrix  of  results 

%  Creating  a  drop  Matrix  using  spherical  coordinates 

[theta,  phi,r]=cart2sph(position_x,position_y,position_z); 

%  Calculating  the  range,  phi,  theta  and  angle  to  every  drop  with 
%  reference  to  phi_0  and  theta_0.  Used  to  estimate  weighted  power  return. 
vector_azimuth=ones(l,drops)*((azimuth*pi)/180); 
vector_elevation=ones(  1  ,drops)*((elevation*pi)/ 180); 
[xf,yf,zfj=sph2cart(vector_azimuth,vector_elevation,r); 
range=sqrt((position_x-xf).A2.+(position_y-yf).A2.+(position_z-zf).A2); 

%  approximation 
angle=atan(range./r); 

%  Saving  the  drop  positions  in  size%d.mat 
%  matrixpos 

matrixpos=[r;phi;theta;angle]; 
save(sprintf('size%d. mat'd), 'matrixpos'); 
end 

% - 

% - 

%  Part  IV,  Calculate  Coherent  Electric  field  return. 

% - 

%  For  every  pulse  the  coherent  Electric  field  return  is  calculatet.  Between 
%  pulses  all  drops  are  moved. 

% - 


%  Run  the  simulation  for  specified  number  of  pulses 

%  Zero  all  vectors  and  matrices  that  are  going  to  be  used 

phase_vector=[]; 

sum_E=[]; 

weighted_E=[]; 

for  p=l  :number_pulses 

disp(sprintf('calculate  the  power-return  of  pulse  %d\n',p)); 
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%  compute  the  (weighted)  I  and  Q  return 

% - 

%  computes  the  total  number  of  drops 

% - 

%  INPUT 
%  -  Matrix  T 
%  -  wavelength 
%  -  angle_3dB 
%  -  distance 
%  -  coord_vol_res 
%  -  number_pulses 

% - 

%  OUTPUT 

%  -  totalE:  the  received  complex  Electric  field) 

%  -  power:  the  power 
%  -  voltage:  the  voltage 

%  -  drops_beam:  the  number  of  drops  in  the  beam 

% - 

%  Only  the  drops  which  are  in  the  resolution  cell  are  selected  for  the 
%  computation 

% - 


number_size=length(matrix_T_Round(l,:)); 

total_E=0; 

total_phase=0; 

new_E=0; 

for  i=T:number_size 

%  opening  the  file  of  registered  data  matrixpos. 
load(sprintf('size%d.maf,i)); 

%  removing  drops  that  are  outside  the  resolution  cell 

number_drops_per_size=length(matrixpos(l,:)); 

dropmatrix=[]; 

%  Filling  the  drop  matrix  with  all  drops  inside  the  beam 

%  resolution  volume 

for  d=l  :number_drops_per_size 

if  ((matrixpos(4,d)<=(angle_3db*pi/180))&(matrixpos(l,d)>=-- 
coord_vol_res(  1 , 1  ))&(matrixpos(  1  ,d)<=coord_vol_res(  1,2))); 
drop  m  a  t  ri  x  ~  [  d  ro  p  m  a  t  r  i  x ,  m  a  t  ri  x  p  o  s  ( : .  d )  [ ; 
end 
end 


% - 

%  Convert  to  cartesian  coordinates 
if  isempty(dropmatrix)==0 

[position_x,position_y,position_z]=sph2cart(dropmatrix(3, :),... 
dropmatrix(2,  :),dropmatrix(  1 , : )); 

%  Create  comparison  vectors 

vector_azimuth=ones(  1  ,length(dropmatrix(  1 ,  :)))*((azimuth*pi)/l  80); 
vector_elevation=ones(  1  dength(dropmatrix(  1  ,:)))*((elevation*pi)/ 180); 
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%  Create  a  phase  vector  for  the  drops  of  current  Diameter 

r_phase=position_x.*u+position_y.*v+position_z.*w; 

phase=exp((2*j*k).*r_phase); 

%  Sum  the  contributions  to  the  E-field 

total_phase=total_phase+sum(phase); 

total_E=total_E+sum(sqrt(matrix_T_Round(2,i)).*phase); 


% - 

%  ceate  weighted  power  return,  concidering  range  (r),  theta,  phi 
%  with  reference  to  boresight  and  radar  parameters 
E_weighted=sqrt(P_t)*G*wavelength/(sqrt(4*pi)A3)./dropmatrix(l,:)... 
.A2.*(exp(-4*log(2).*((dropmatrix(3,:)-vector_azimuth).A2/... 
(2*angle_3db*pi/180)A2+((dropmatrix(2,:)-vector_elevation)... 
,A2/(2*angle_3db*pi/180)A2)))); 

new_E=new_E+sum(E_weighted.*(sqrt(matrix_T_Round(2,i)).*phase)); 

end 

end 

% - 


% - 

%  This  part  simulates  the  movement  of  the  drops  and  gives  new  positions. 

% - 

%  INPUT 

%  -  matrix  T  Round 
%  -  PRT 
%  -  azimuth 
%  -  elevation 
%  -  rain_anglel 
%  -  rain_angle2 
%  -  diam_lim 
%  -  resolution  1 
%  -  resolution2 

%  -  wind_vector 

% - 

%  OUTPUTS 

%  Output  files  containing  sizes  with  new  positions 
%  structured  the  same  way  as  in  previous 

%  there  is  also  an  output  file  for  speed  (one  for  each  size)  for  the  pulse 
%  pair. 

% - 


number_size=length(matrix_T_Round(l,:)); 

%  adding  a  Gaussian  random  term  to  make  the  wind  speed  spread 
if  wind_vector(l)~=0 

new_wind_vector(  1  )=wind_vector(  1  )*normrnd(  1  ,abs(spread/max(wind_vector(  1 )))); 
else 

new_wind_vector(  1  )=0; 
end 

if  wind_vector(2)~=0 

new_wind_vector(2)=wind_vector(2)*normrnd(l,abs(spread/max(wind_vector(2)))); 

else 

new_wind_vector(2)=0; 

end 
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if  wind_vector(3)~=0 

new_wind_vector(3)=wind_vector(3)*normrnd(l,abs(spread/max(wind_vector(3)))); 

else 

new_wind_vector(3)=0; 

end 

for  i=l  :number_size 

%  opening  of  data  file  registered  in  matrixpos. 
load(sprintf('size%d.maf,i)); 

%  total  number  of  drops 

drops=matrix_T_Round(3,i); 

%  transforming  to  rectangular  coordinates 
[x,y,z]=sph2cart(matrixpos(3,:),matrixpos(2,:),matrixpos(1 ,:)); 

diameter=matrix_T_Round(  1  ,i); 


% - 

%  Drop  fall  speed  depends  on  drop  diameter  in  the  z  axis, 

%  Atlas-Ulbrich  approximation  is  used 

%  Approximation  is  valid  in  the  diameter  range  5*e-4,  5*e-3 
% - 


wtmax=386.6*((matrix_T_Round(l,i)*0.001)A0.67); 

velocityz=(wtmax*(ones(l,drops)))*cos((rain_anglel*pi)/180); 

velocityx=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 

cos((rain_angle2*pi)/180); 

velocityy=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l  80); 


newest_wind_vector=[] ; 


%  adding  Gaussian  random  term  to  make  the  wind  speed 
%  spread  between  drops 
if  new_wind_vector(l)~=0 

ne  west_wind_vector(  1 , :  )=new_wind_vector(  1 )  * . . . 

normrnd(l,abs(spread_2/max(new_wind_vector(l))),l, drops); 


else 

newest_wind_vector(l,:)=0; 

end 

if  new_wind_vector(2)~=0 

newest_wind_vector(2,:)=new_wind_vector(2)*... 

normrnd(l,abs(spread_2/max(new_wind_vector(2))),l, drops); 


else 

newest_wind_vector(2,:)=0; 

end 

if  new_wind_vector(3)~=0 
ne  west_wind_vector(3 , :  )=new_wind_vector(  3 )  * . . . 

normrnd(l,abs(spread_2/max(new_wind_vector(3 ))),!, drops); 


else 

ne  west_wind_vector(  3 , :  )=0; 
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end 


%  If  spread  within  pulse 
velocityxt=newest_wind_vector(l,:)-velocityx; 
velocityyt=newest_wind_vector(2,:)-velocityy; 
velocityzt=newest_wind_vector(3,:)-velocityz; 


%  speed  registration 

m_vclocity— |  velocityxt;vclocityyt;velocityzt  |; 
save(sprintf('velocity%d.mat',i),'m_velocity'); 

deplacementx=velocityxt*PRT ; 
deplacementy=velocityyt*PRT ; 
deplacementz=velocityzt*PRT ; 

position_x=x+deplacementx; 

position_y=y+deplacementy; 

position_z=z+deplacementz; 

%  transformation  into  spherical  coordinates 

[theta,  phi,r]=cart2sph(position_x,position_y,position_z); 

%  computation  of  new  angles 
vector_azimuth=ones(l,drops)*(azimuth*pi)/180; 
vector_elevation=ones(l,drops)*(elevation*pi)/180; 
[xf,yf,zf]=sph2cart(vector_azimuth,vector_elevation,r); 

range=sqrt((position_x-xf).A2.+(position_y-yf).A2.+(position_z-zf).A2); 
%  approximation 
angle=atan(range ./ r) ; 

%  matrixpos 

matrixpos=[r;phi;theta;angle]; 

%  save  to  file,  matrixpos 

s=sprintf(  'size%d.  mat' ,  i) ; 
save(s,  'matrixpos'); 


end 

%■ 


%  add  to  E-field  and  phase  vector 
sum_E=[sum_E,total_E] ; 
weighted_E=[weighted_E,new_E]; 
phase_vector=[phase_vector,ANGLE(total_phase)]; 

end 


%• 


% - 

%  Power  return  for  every  pulse  including  error 

running_Z=(Norm*abs(sum_E).A2)*wavelengthA4/piA5/abs(kw)A2/volume_resolution... 

/le-18+error_Z; 
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power_retum=Norm*abs(weighted_E).A2; 


%  [W] 


%  calculate  average  Power  and  Z 
Z_avg=sum(running_Z)/p; 

Z_avg_plot=ones(  1  ,p)*Z_avg; 

power_return_avg=sum(Norm*abs(weighted_E).A2)/p;  %  [W] 

power_return_avg_plot=ones(l,p)*power_return_avg; 

%  calculate  doppler  frequency 
phase_vector_u=unwrap(phase_vector); 

doppler_fq=([phase_vector_u(2:p)-phase_vector_u(l  :p- 1  )])/2/pi/PRT ; 
number_vector=l  :p; 


% - 

%  Calculate  references 
%  Integral  form  of  Z 
ii=0: 100/1000: 100; 

zz=ii.A6.*8000.*exp(-4.1.*rain_rate.A(-.21).*ii); 

real_Z=ones(l,number_pulses)*10*logl0(trapz(ii,zz)); 

%  reflectivity_real=ones(l,number_pulses)*10*logl0(piA5/wavelengthA4*... 

%  volume_resolution*abs(kw)A2*real_Z*le-18); 

%  Adding  the  error  to  the  weighted_E  estimations 

error_power=error_rcs*P_t*GA2*wavelengthA2/piA2*(2*angle_3db*pi/180)A2*... 

(c*pulsewidth*le-6)/1024/log(2)/piA2/(distance*le3)A2/volume_resolution; 

%  Evaluate  average  Power  return  to  make  Z  estimation  including  error 
Z_estimation=power_return_avg/P_t/GA2/wavelengthA2*... 

piA2/(2*angle_3db*pi/180)A2/(c*pulsewidth*le-6)*1024*log(2)*(distance*le3)A2/... 

((kwA2)*(piA5)/(wavelengthA4))/le-18; 

Z_estimation_plot=ones(l,number_pulses)*10*logl0(Z_estimation+error_Z); 

%  calculate  the  FFT  of  the  weighted  electric  field  return  (1024  points)  to 
%  get  the  power  spectal  estimate 
NN=1024;  %  number  of  FFT  points 

spectral_estimate=abs(fft(weighted_E-mean(weighted_E),NN)).A2; 

spectral_estimate_shifted=fftshift(spectral_estimate); 

% - 

%  Part  V,  Plot  of  results 

% - 


figure)  1); 

subplot(2,l,l); 

plot(number_vector(l  :p-l),doppler_fq(l:p-l),'-x',number_vector(  1  :p),... 
ones(l,p)*mean(doppler_fq),'— r'); 

title(['Doppler  frequency  for  all  pulse  pairs  for  PRF=  ',num2str(l/PRT),'  Hz']); 

xlabel('Pulse  pair'); 

ylabel('fd'); 

legend('Doppler  frequency  per  pulse  pair'. 'Average  Doppler  estimation') 
grid  on 


subplot(2,l,2); 

v_r=doppler_fq*wavelength/2 ; 
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hist(v_r,  15) 
hold  on 

[ff,xx]  =  ksdensity(v_r); 

plot(xx,ff/max(ff)*max(hist(v_r,15))) 

title('Radial  Velocity  Histogram  with  Gaussian  Fit'); 

xlabel('Radial  velocity  [m/s]'); 

ylabel('Frequency'); 

grid  on 

figure  (2) 

subplot(2,l,l) 

%  also  converting  to  dBm 

plot(number_vector,10*logl0(power_return+error_power)+30,'-x',... 

number_vector,10*logl0(power_return_avg_plot+error_power)+30,'— r'); 
title('Power  Return  for  all  pulses'); 
xlabel('Sample  time'); 
ylabel('Power  [dBm]'); 

legend('Power  return',  'Average  Power  return', 0) 
grid  on 

subplot(2,l,2); 

plot(number_vector*1000*PRT,10*logl0(abs(running_Z)),'-xb',number_vector... 
*1000*PRT,10*logl0(abs(Z_avg_plot)),'-.r',number_vector*1000*... 
PRT,real_Z,'-k',number_vector*1000*PRT,Z_estimation_plot,':b'); 
title('Estimated  Z  per  pulse'); 
xlabel('Sample  time'); 
ylabel('Z  [dBZ]'); 

legend('Simulation  Z,  actual',  'Average  Z','Z=\int  DA6N(D)dD',... 

'Estimated  Z,  Zeroth  moment', 0) 
hline  =  findobj(gca,'Type','line'); 
set(hline,'LineWidth',2) 
grid  on 


figure  (3) 
subplot(2,l,l) 

f=abs(xcorr((weighted_E),'coeff)); 

tau_ts=ones(l,length(f))*f((length(f)+l)/2-l); 

plot(l:length(f),abs(xcorr((weighted_E),'coeff)),'b',l:length(f),tau_ts,'— r'); 
title('Return  Voltage  Correlation'); 
xlabel('Number  of  pulses'); 
ylabel('Normalized  correlation  value'); 

legend(['Correlation  function'], [Vho(\tau=T_s)=',num2str(tau_ts(l))]) 
grid  on 

subplot(2,l,2) 

%  Prepare  for  ticks  so  every  100  Hz  is  marked 
n=0:l/PRT/100; 

N_i=100*(NN-l)*n*PRT; 
alfa=zeros(  1 , 1  /PRT/ 1 00+ 1 ); 
vic=(length(alfa)- 1  )/2; 
for  beta=l  dength(alfa) 
alfa(  1  ,beta)=-vic*  100; 
vic=vic-l; 
end 
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a  \'g_t  q  =mean(d  o  p  p  I  c  r_fq ) 
std_fq=std(dopp  1  erfq) 
gauss_vectoi -min(  alfa) :  1 0 :  max(  alfa) ; 
uu=length(gauss_vector); 

[value  1,  index  l]=max(abs(spectral_estimate_shifted)); 
avg_fq_fft=(  index  1  -NN/2)/PRT /NN ; 

plot_gauss_vector=linspace(0,NN,uu); 

plot_norm_fq=normpdf(gauss_vector,avg_fq,std_fq/2); 

plot_norm_fft=normpdf(gauss_vector,avg_fq_fft,std_fq/2); 

plot(spectral_estimate_shifted/max(spectral_estimate_shifted)) 
hold  on 

plot(fftshift(abs(fft(xcorr((weighted_E),'coeff),1024))/... 

max(abs(fft(xcorr((weighted_E),'coeff ),  1 024)))), '-k') 
title('Power  Spectral  Estimate'); 
xlabel('Doppler  Frequancy  [Hz]'); 
ylabel('S|f|,  normalized'); 
xlim([0  NN]); 
set(gca,'xtick',  N_i) 
set(gca,'xticklabel',  num2str(alfa')); 
grid  on 

2)  WR_dual_prf 

%  Weather  Radar  Signal  Simulator,  developed  version 
%  Two  PRFs  enabling  pulse  pair  comparison 
%  Ulf  Schroder 
%  August,  2005 

% - 


%  Main  program  for  The  NPS/Marielle  Gosset  Weather  Radar  Signal  Simulator. 


% - 

%  DESCRIPTION 
%  Part  I 

%  From  Radar  Data: 

%  -  Generate  Beam  Resolution  cell  approximation 

%  -  Generate  Cube  adding  a  margine  to  the  Beam  Resolution  cell  to  allow  for 
%  the  rain  drops  to  fall  trough  the  through  beam  over  time. 

% 

%  Part  II 

%  From  Rain  Parameters: 

%  Generating  the  fundamental  T  Matrix  concisting  of 
%  -  Diameter  Vector  and  Number  of  drops  in  each  diameter  interval 
%  using  the  Marshall-Palmer  exponential  approximation. 

%  -  the  backscatter  RCSs  for  each  Diameter  size 
%  -  the  total  number  of  drops  for  each  size  per  mA3 

%  Part  III 

%  Randomly  place  all  drops  in  cube 
%  Part  IV 

%  For  every  pulse  the  coherent  Electric  field  return  is  calculated.  Between 
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%  pulses  all  drops  are  moved. 


%  Part  V 

%  Calculate  and  plot  results 

% - 

clear  all 
close  all 
clc 

%  Input  Data 

c=3e8;  %  Speed  of  light  [m/s] 

%  Radar  Data 

P_t=  1 0e3 ;  %  T ransmitted  power 

G_dB=20;  %  Maximum  gain  in  dB 

G=10A(G_dB/10);  %  Maximum  gain 

angle_3db=0.20;  %  3  dB  Beam  angle  [degrees],  (half  half  power  beamwidth) 

pulsewidth=0.20;  %  Pulsewidth  [microsecond] 

elevation=3;  %  Elevation  angle  [degrees] 

margine=5;  %  margin  for  resolution  cell[m] 

azimuth=0;  %  Azimuth  angle  [rad] 

number_pulses=40;  %  Number  of  pulses 

PRT_l=l/2000;  %  Doppler  Pulse  Repetition  Time  [seconds] 

PRT_2=l/200;  %  Avg  Power  Pulse  Repetition  Time  [seconds] 

distance=0.5;  %  Range  from  the  radar  to  the  resolution  cell  [Km] 
wavelength=0.1;  %  Wavelength  [m] 

k=2*pi/wavelength; 

%  Rain  parameters 

diam_low=0.05;  %  Lowest  diameter  of  Precipitation  [mm] 

diam_lim=l;  %  Divider  between  low  and  mid  [mm] 

diam_mid=5;  %  Divider  between  mid  and  high  [mm] 

diam_high=15;  %  Highest  diameter  of  Precipitation  [mm] 

resolutionl=0.08;  %  Lower  diameter  resolution  [mm] 
resolution2=0.3;  %  Higher  diameter  resolution  [mm] 
resolution3=0.8;  %  Highest  diameter  resolution  [mm] 
sort_approx=l;  %  RCS  approximation  model;  [l]=Rayleigh 
kw=0.93;  %  abs(K)A2  approx  for  water 

rain_rate=50;  %  Rain  Rate  [mm/h] 

rain_anglel=0;  %  degrees  Rain  Angle  (Phi) 

rain_angle2=30;  %  degrees  Rain  Angle  (Theta) 
wind_vcctoi'-| - 1  ();();()[;  %  Wind  speed  vector  (x,y,z)  [m/s] 
spread=l;  %  Wind  speed  spread  [m/s],  (all  directions) 

spread_2=4;  %  Drop  velocity  spread  [m/s], 

Normal 00;  %  Normalization  coefficient 

%  Directional  cosines 

u=sin(pi/2-elevation*pi/180)*cos(azimuth*pi/180); 

v=sin(pi/2-elevation*pi/180)*sin(azimuth*pi/180); 

w=cos(pi/2-elevation*pi/180); 

% - 

%  Part  I,  Generate  the  resolution  volume 

% - 

%  Approximation  of  the  resolution  cell  of  the  pulse  circular  box, 
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%  for  any  azimuth  angle 

% - 

%  INPUT 
%  -  distance 
%  -  angle_3dB 
%  -  pulsewidth 
%  -  elevation 
%  -  margine 
%  -  azimuth 

% - 

%  OUTPUT 

%  -  volume_resolution:  computed  volume  of  the  resolution  cell 
%  -  volume_box:  volume  of  the  box  where  the  drops  are  put 
%  -  box:  the  8  rectangular  coordinates  of  the  box 
%  -  coord_vol_res:  the  coordinates  of  the  resolution  cell  where  the  drops 
%  are  taken. 

% - 


% - 

%  Creating  the  circular  beam  resolution  cell 
range_res=(c*pulsewidth*le-6)/2; 

%  Range  resolution  (ct/2)  [m] 

%  small  angle  approximation  allows  for  following 
angle_rad=(angle_3  db  *pi)/ 1 8  0 ; 

%  3  dB  beamwidth  [rad] 
radius  l=(tan(angle_rad)*(distance*  1000)); 

%  Radius  of  beam  at  R=distance 
radius2=(tan(angle_rad)*((distance*  1000)+range_res));  ... 

%  Radius  of  beam  at  R=distance  +  range  resolution 

surfacel=pi*radiuslA2; 

surface2=pi*radius2A2; 

%  volume  of  cone  =  base  surface  times  height  for  all  three 
volume_conel=(surfacel*(distance*1000))/3; 
volume_cone2=(surface2*((distance*1000)+range_res))/3; 
volume_resolution=(volume_cone2-volume_conel);  %  Beam  Cone  Volume 

%  computation  of  the  resolution  cell  coordinates,  spherical  coordinates 
%  elevation  in  radians 

elevation_rad=(elevation*pi)/180;  %  Elevation  [rad] 

azimuth_rad=(azimuth*pi)/180;  %  Azimuth  [rad] 

%  Matrix  used  when  comparing  wheather  a  drop  is  inside  or  outside 
%  beamwidth 

coord_vol_res=[distance*  1000, (distance*  1000)+range_res;elevation_rad-angle_rad,... 

elevation_rad+angle_rad;azimuth_rad-angle_rad,azimuth_rad+angle_rad]; 

% - 


% - 

%  Creating  the  box 

%  computation  of  the  box  and  its  coordinates 

coin  1 =[azimuth_rad+angle_rad,elevation_rad+angle_rad,  distance*  1 000] ; 
coin2~[  azimuth_rad-angle_rad.clevation_rad^angle_rad. distance*  1 000 1; 
coin3=[azimuth_rad+angle_rad,elevation_rad-angle_rad,distance*1000]; 
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coin4=[azimiith_rad-angle_rad,elevation_rad-angle_rad,distance*1000]; 
coin5=[azimuth_rad+angle_rad,elevation_rad+angle_rad,distance*1000+range_res]; 
coin6=[azimuth_rad-angle_rad, elevation_rad+angle_rad, distance*  1000+range_res]; 
coin7=[azimuth_rad+angle_rad,elevation_rad-angle_rad, distance*  1000+range_res]; 
coin8=[azimuth_rad-angle_rad,elevation_rad-angle_rad, distance*  1000+range_res]; 

[x  1  ,y  1  ,zl  ]=sph2cart(coin  1  ( 1  ),coin  1  (2), coin  1  (3 )); 
[x2,y2,z2]=sph2cart(coin2(l),coin2(2),coin2(3)); 

[x3  ,y3  ,z3]=sph2cart(coin3  ( 1  ),coin3  (2),coin3  (3 )); 
[x4,y4,z4]=sph2cart(coin4(l),coin4(2),coin4(3)); 
[x5,y5,z5]=sph2cart(coin5(l),coin5(2),coin5(3)); 
[x6,y6,z6]=sph2cart(coin6(l),coin6(2),coin6(3)); 
[x7,y7,z7]=sph2cart(coin7(l),coin7(2),coin7(3)); 
[x8,y8,z8]=sph2cart(coin8(l),coin8(2),coin8(3)); 

xmax=max([xl,x2,x3,x4,x5,x6,x7,x8]); 
xmin=min([xl,x2,x3,x4,x5,x6,x7,x8]); 
ymax=max([y  1  ,y2,y3,y4,y5,y6,y7  ,y8]); 
ymin=min([y  1  ,y2,y3,y4,y5,y6,y7  ,y8]); 
zmax=max([zl,z2,z3,z4,z5,z6,z7,z8]); 
zmin=min([zl,z2,z3,z4,z5,z6,z7,z8]); 

%  the  box  with  its  margin 

coin_l=[xmin-margine,ymin-margine,zmin-margine]; 

coin_2=[xmin-margine,ymax+margine,zmin-margine]; 

coin_3=[xmax+margine,ymax+margine,zmin-margine]; 

coin_4=[xmax+margine,ymin-margine,zmin-margine]; 

coin_5=[xmin-margine,ymin-margine,zmax+margine]; 

coin_6=[xmin-margine,ymax+margine,zmax+margine]; 

coin_7=[xmax+margine,ymax+margine,zmax+margine]; 

coin_8=[xmax+margine,ymin-margine,zmax+margine]; 


box=[coin_r,coin_2',coin_3',coin_4',coin_5',coin_6',coin_7',coin_8']; 

%  calcidate  the  volume  of  the  box 

volume_box=abs((xmax+margine-(xmin-margine))*(ymax+margine-(ymin-margine))*(... 

zmax+margine-(zmin-margine))); 


cube=[coin_r,coin_2',coin_3',coin_4',coin_r,coin_5',coin_6',coin_2',... 

coin_6',coin_7',coin_3',coin_7',coin_8',coin_4',coin_8',coin_5']; 

% - 


% - 

%  Part  II,  Generate  the  fundamental  T  Matrix 

% - 

%  This  part  produces  a  matrix  of  all  drops  that  contribute  to  the  total 
%  reflectivity  of  the  resolution  cell 

% - 

%  First  all  sizes  are  created,  with  two  vectors  of  different  size  resolution. 

%  The  size  vectors  represents  all  diameters.  Next  the  backscatter  RCS  for  each 
%  size  is  computed. 

%  Using  the  Marshall-Palmer  exponential  approximation,  the  number  of  drops  of 
%  each  diameter  is  computed. 

% - 
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%  INPUT 
%  -  diam_low 
%  -  diam_lim 
%  -  diam_high 
%  -  resolution  1 
%  -  resolution2 
%  -  sort_approx 
%  -  wavelength 
%  -  Kw 
%  -  rain_rate 

% - 

%  OUTPUT 
%  -  number_drops 
%  Matrix  T : 

%  -  first  row  all  sizes 
%  -  second  row  all  the  backscatter  RCSs 
%  -  third  row  the  total  number  of  drops  for  each  size  per  mA3 
%  the  columns  start  with  the  smallest  size 

% - 

%  Create  dropsize  vector 

%  Precipitaion  Diameter  vector  (lowest<D<divider) 
vector_size_l=diam_low:resolutionl:diam_lim; 

%  Precipitaion  Diameter  vector  (divider<D<mid) 
vector_size_2=diam_lim:resolution2:diam_mid; 

%  Precipitaion  Diameter  vector  (mid<D<highest) 
vector_size_3=diam_mid:resolution3:diam_high; 

%  Precipitaion  Diameter  Vector 

vector_size=[vector_size_l,vector_size_2,vector_size_3]; 


% - 

%  Create  RCS  vector  for  all  Rain  drops  of  Precipitation  Vector. 
rcs=((kwA2)*(piA5)/(wavelengthA4))*((vector_size*le-3).A6); 
% - 


% - 

%  The  Marshall  Palmer  approximation  is  used  to  create  Precipitaion  Diameter 
%  Distribution  Vector  (Drop  Size  Distribution  N(D)  Vector) 

%  The  number  of  drops  are  scaled  down  by  a  factor  (Norm)  to  improve  the 
%  speed  of  the  simulator 

% - 

lambda=4.1*(rain_rateA(-0.21));  %[mm-l] 

number=8000*exp(-lambda*vector_size)/Norm; 

vector_resl=ones(l,length(vector_size_l))*resolutionl; 

vector_res2=ones(l,length(vector_size_2))*resolution2; 

vector_res3=ones(l,length(vector_size_3))*resolution3; 

%  delta(Diameter)  dD 

vector_res=[vector_resl,vector_res2,vector_res3]; 

%  the  approximate  "true"  number  of  drops  N(D)dD  vector 
number=number.*vector_res; 


99 


%• 


% - 

%  create  matrix  T 
T=[vector_size;rcs;number]; 


% - 

%  Change  matrix  T  to  matrix  'matrixTRound' 

% - 

%  This  part  computes  the  total  number  of  drops  by  utilizing  the  resolution 
%  cell  and  the  matrix_T_Round 

%  The  only  difference  between  the  matrix  T  Round  and  T,  is  the  fact  that 
%  matrix_T_Round  has  the  total  number  of  drops  on  the  third  row  and  they 
%  are  rounded  off  to  nearest  integer  value. 

% - 


matrix_T_Round=T ; 

matrix_T_Round(3,:)=round(matrix_T_Round(3,:)*volume_box); 

%  Keep  track  of  errors  due  to  round  off 
diff=T(3 , :  )*volume_box-matrix_T_Round(3, : ); 

error_rcs=sum(diff.*matrix_T_Round(2,:))*Norm/volume_box*volume_resolution; 

error_Z=error_rcs/volume_resolution*wavelengthA4/piA5/abs(kw)A2/le-18; 


% - 

%  Part  III,  Initial  positioning  of  the  drops  within  the  volume 

% - 

%  This  part  generates  the  initial  position  of  the  drops  in  the  volume  of 
%  the  box. 

%  The  positions  are  recorded  in  the  file  "sizexx.maf  taking  into  account 
%  the  azimuth  angle 

% - 

%  INPUT 
%  -  box 

%  -  matrix_T_Round 
%  -  elevation 
%  -  azimuth 

% - 

%  OUTPUT 

%  All  files  including  the  sizes.  Each  contains  the  positions  of  all  drops 
%  for  a  given  diameter.  The  data  recoreded  as  follows: 

%  -  radial  distance  r 
%  -  elevation  phi 
%  -  azimuth  theta 

%  -  angle  with  respect  to  beam  center.  This  parameter  will  be 
%  used  to  identify  drops  which  are  in  the  beam. 

% - 


%  Randomly  put  the  calculated  number  of  drops  of  each  diameter  in  the  box 
number_size=length(matrix_T_Round(l,:)); 

for  i=l  :number_size 

%  select  number  of  drops  of  certain  Diameter 
drop  s=matrix_T_Round(  3 ,  i ) ; 
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%  create  as  many  positions  on  the  x  axis  as  #  drops  per  size 
elements_x=rand(  1 ,  drops); 

%  put  them  in  the  range  of  interest 
%  Random  placement  x  for  every  drop 

position_x=box(  1 , 1  )+(box(  1 ,3  )-box(  1 , 1  ))*elements_x; 

%  create  as  many  positionson  the  y  axis  as  #  drops  per  size 
elements_y=rand(  1  ,drops); 

%  put  them  in  the  range  of  interest 
%  Random  placement  y  for  every  drop 

position_y=box(2, 1  )+(box(2,2)-box(2, 1  ))*elements_y; 

%  create  as  many  positions  on  the  z  axis  as  #  drops  per  size 
elements_z=rand(  1 , drops); 

%  put  them  in  the  range  of  interest 
%  Random  placement  z  for  every  drop 

position_z=box(3,l)+(box(3,5)-box(3,l))*elements_z; 

%  matrix  of  results 

%  Creating  a  drop  Matrix  using  spherical  coordinates 

[theta,  phi,r]=cart2sph(position_x,position_y,position_z); 

%  Calculating  the  range,  phi,  theta  and  angle  to  every  drop  with 
%  reference  to  phi_0  and  theta_0.  Used  to  estimate  weighted  power  return. 
vector_azimuth=ones(l,drops)*((azimuth*pi)/180); 
vector_elevation=ones(  1  ,drops)*((elevation*pi)/180); 
[xf,yf,zf]=sph2cart(vector_azimuth,vector_elevation,r); 
range=sqrt((position_x-xf).A2.+(position_y-yf).A2.+(position_z-zf).A2); 

%  approximation 
angle=atan(range./r); 

%  Saving  the  drop  positions  in  size%d.mat 
%  matrixpos 

matrixpos=[r;phi;theta;angle]; 
save(sprintf('size%d. mat'd), 'matrixpos'); 
end 

% - 


% - 

%  Part  IV,  Calculate  Coherent  Electric  field  return. 
% - 


%  For  every  pulse  the  coherent  Electric  field  return  is  calculatet.  Between 
%  pulses  all  drops  are  moved. 


%■ 


%  Run  the  simulation  for  specified  number  of  pulses 

%  Zero  all  vectors  and  matrices  that  are  going  to  be  used 

powerTM); 

vector_power=[]; 
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phase_vector=zeros(  1 ,2*number_pulses); 
imaginary_retiirn=zeros(  1 ,2*number_pulses); 
sum_E=[]; 
weighted_E=[]; 

for  p=l  :number_pulses 

disp(sprintf('calculate  the  power-return  of  pulse  %d\n',p)); 

%  compute  the  (weighted)  I  and  Q  return 

% - 

%  computes  the  total  number  of  drops 

% - 

%  INPUT 
%  -  Matrix  T 
%  -  wavelength 
%  -  angle_3dB 
%  -  distance 
%  -  coord_vol_res 
%  -  number_pulses 

% - 

%  OUTPUT 

%  -  total  E:  the  received  complex  Electric  field) 

%  -  power:  the  power 
%  -  voltage:  the  voltage 

%  -  drops_beam:  the  number  of  drops  in  the  beam 

% - 

%  Only  the  drops  which  are  in  the  resolution  cell  are  selected  for  the 
%  computation 

% - 


number_size=length(matrix_T_Round(l,:)); 

total_E=0; 

total_phase=0; 

new_E=0; 

for  i=l:number_size 

%  opening  the  file  of  registered  data  matrixpos. 
load(sprintf('size%d. mat',i)); 

%  removing  drops  that  are  outside  the  resolution  cell 

number_drops_per_size=length(matrixpos(l,:)); 

dropmatrix=[]; 

%  Filling  the  drop  matrix  with  all  drops  inside  the  beam 

%  resolution  volume 

for  d=l  :number_drops_per_size 

if  ((matrixpos(4,d)<=(angle_3db*pi/180))&(matrixpos(l,d)>=-- 
coord_vol_res(  1 , 1  ))&(matrixpos(  1  ,d)<=coord_vol_res(  1,2))); 
dropmatrix=[dropmatrix,matrixpos(:,d)]; 
end 
end 


% - 

%  Convert  to  cartesian  coordinates 
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if  isempty(dropmatrix)==0 

[position_x,position_y,position_z]=sph2cart(dropmatrix(3, 
dropmatrix(2,  :),dropmatrix(  1 , : 

%  Create  comparison  vectors 

vector_azimuth=ones(l,length(dropmatrix(l,:)))*((azimuth*pi)/180); 
vector_elevation=ones(  1  dength(dropmatrix(  1  ,:)))*((elevation*pi)/ 180); 

%  Create  a  phase  vector  for  the  drops  of  current  Diameter 

r_phase=position_x.*u+position_y.*v+position_z.*w; 

phase=exp((2*j*k).*r_phase); 

%  Sum  the  contributions  to  the  E-field 

total_phase=total_phase+sum(phase); 

total_E=total_E+sum(sqrt(matrix_T_Round(2,i)).*phase); 


% - 

%  ceate  weighted  power  return,  concidering  range  (r),  theta,  phi 
%  with  reference  to  boresight  and  radar  parameters 

E_weighted=sqrt(P_t)*G*wavelength/(sqrt(4*pi)A3)./dropmatrix(l,:).A2.*... 
(exp(-4*log(2).*((dropmatrix(3,:)-vector_azimuth).A2/... 
(2*angle_3db*pi/180)A2+((dropmatrix(2,:)-vector_elevation)  A2/... 
(2*angle_3db*pi/180)A2)))); 

new_E=new_E+sum(E_weighted.*(sqrt(matrix_T_Round(2,i)).*phase)); 

end 

end 

phase_vector(2*p- 1  )=ANGLE(total_phase); 
imaginary_return(2*p-l)=total_phase; 


yo¬ 


yo - 

%  This  part  simulates  the  movement  of  the  drops  and  gives  new  positions. 
%  First  using  the  Doppler  PRT 

% - 

%  INPUT 

%  -  matrix  T  Round 
%  -  PRT 
%  -  azimuth 
%  -  elevation 
%  -  rain_anglel 
%  -  rain_angle2 
%  -  diam_lim 
%  -  resolution  1 
%  -  resolution2 

%  -  wind_vector 

% - 

%  OUTPUTS 

%  Output  files  containing  sizes  with  new  positions 
%  structured  the  same  way  as  in  previous 

%  there  is  also  an  output  file  for  speed  (one  for  each  size)  for  the  pulse 
%  pair. 

% - 
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number_size=length(matrix_T_Round(l,:)); 

%  adding  a  Gaussian  random  term  to  make  the  wind  speed  spread 
if  wind_vector(l)~=0 

new_wind_vector(  1  )=wind_vector(  l)*normrnd(  1  ,abs(spread/max(wind_vector(  1 )))); 
else 

new_wind_vector(  1  )=0; 
end 

if  wind_vector(2)~=0 

new_wind_vector(2)=wind_vector(2)*normrnd(l,abs(spread/max(wind_vector(2)))); 

else 

new_wind_vector(2)=0; 

end 

if  wind_vector(3)~=0 

new_wind_vector(3)=wind_vector(3)*normrnd(l,abs(spread/max(wind_vector(3)))); 

else 

new_wind_vector(3)=0; 

end 


for  i=l  :number_size 

%  opening  of  data  fde  registered  in  matrixpos. 
load(sprintf('size%d.mat',i)); 

%  total  number  of  drops 
drops=matrix_T_Round(3,i); 

%  transforming  to  rectangular  coordinates 
[x,y,z]=sph2cart(matrixpos(3,:),matrixpos(2,:),matrixpos(1 

diameter=matrix_T_Round(  1  ,i); 


% - 

%  Drop  fall  speed  depends  on  drop  diameter  in  the  z  axis, 

%  Atlas-Ulbrich  approximation  is  used 

%  Approximation  is  valid  in  the  diameter  range  5*e-4,  5*e-3 
% - 


wtmax=386.6*((matrix_T_Round(l,i)*0.001)A0.67); 

velocityz=(wtmax*(ones(l,drops)))*cos((rain_anglel*pi)/180); 

velocityx=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 

cos((rain_angle2*pi)/180); 

velocityy=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l  80); 

newest_wind_vector=[]; 


%  adding  Gaussian  random  term  to  make  the  wind  speed 
%  spread  between  different  drops  sizes 
if  new_wind_vector(l)~=0 

ne  west_wind_vector(  1 , :  )=new_wind_vector(  1 )  *normrnd(  1 , . . . 
abs(spread_2/max(new_wind_vector(l))),l, drops); 

else 
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newest_wind_veetor(l,:)=0; 

end 

if  new_wind_vector(2)~=0 

newest_wind_vector(2,:)=new_wind_vector(2)*normrnd(l,... 
abs(spread_2/max(new_wind_vector(2))),l,  drops); 

else 

newest_wind_vector(2,:)=0; 

end 

if  new_wind_vector(3)~=0 

newest_wind_vector(3,:)=new_wind_vector(3)*normrnd(l,... 
abs(spread_2/max(new_wind_vector(3))),l, drops); 

else 

ne  west_wind_vector(  3 , :  )=0; 
end 

velocityxt=newest_wind_vector(  1 , :  )-velocityx; 
velocityyt=newest_wind_vector(2,:)-velocityy; 
velocityzt=newest_wind_vector(3,:)-velocityz; 

%  speed  registration 

m_velocity=[velocityxt;velocityyt;velocityzt]; 

save(sprintf('velocity%d.maf,i),'m_velocity'); 

deplacementx=velocityxt*PRT_l ; 
deplacementy=velocityyt*PRT_l ; 
deplacementz=velocityzt*PRT_l ; 

position_x=x+deplacementx; 

position_y=y+deplacementy; 

position_z=z+deplacementz; 

%  transformation  into  spherical  coordinates 

[theta,  phi,r]=cart2sph(position_x,position_y,position_z); 

%  computation  of  new  angles 
vector_azimuth=ones(l,drops)*(azimuth*pi)/180; 
vector_elevation=ones(l,drops)*(elevation*pi)/180; 
[xf,yf,zf]=sph2cart(vector_azimuth,vector_elevation,r); 

range=sqrt((position_x-xf).A2.+(position_y-yf).A2.+(position_z-zf).A2); 

%  approximation 

angle=atan(range./r); 

%  matrixpos 

matrixpos=[r;phi;theta;angle]; 

%  save  to  fde,  matrixpos 

s=sprintf(  'size%d.  mat' ,  i) ; 
save(s,  'matrixpos'); 


end 


% - 

%  Calculate  phase  for  pulse  pair 
total_phase=0; 
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for  i=l: number  size 


%  opening  the  file  of  registered  data  matrixpos. 
load(sprintf('size%d. mat',i)); 

%  removing  drops  that  are  outside  the  resolution  cell 
number_drops_per_size=length(matrixpos(  1 ,:)); 
dropmatrix=[]; 

%  Filling  the  drop  matrix  with  all  drops  inside  the  beam 

%  resolution  volume 

for  d=l  :number_drops_per_size 

if  ((matrixpos(4,d)<=(angle_3db*pi/180))&(matrixpos(l,d)>=... 
coord_vol_res(  1 , 1  ))&(matrixpos(  1  ,d)<=coord_vol_res(  1,2))); 
dropmatrix=[dropmatrix,matrixpos(:,d)]; 
end 
end 


% - 

%  Convert  to  cartesian  coordinates 
if  isempty(dropmatrix)==0 

[position_x,position_y,position_z]=sph2cart(dropmatrix(3, :),... 
dropmatrix(2,  :),dropmatrix(  1 , : )); 

%  Create  comparison  vectors 

vector_azimuth=ones(ldength(dropmatrix(l,:)))*((azimuth*pi)/180); 
vector_elevation=ones(  1  Jength(dropmatrix(  1  ,:)))*((elevation*pi)/ 180); 


%  Create  a  phase  vector  for  the  drops  of  current  Diameter 

r_phase=position_x.*u+position_y.*v+position_z.*w; 

phase=exp((2*j*k).*r_phase); 

%  Sum  the  contributions  to  the  E-field 
total_phase=total_phase+sum(phase); 
end 
end 

phase_vector(2*p)=ANGLE(total_phase); 

imaginary_return(2*p)=total_phase; 


% - 

%  This  part  simulates  the  movement  of  the  drops  and  gives  new  positions 
%  using  the  Avg  Power  PRT 

% - 


number_size=length(matrix_T_Round(l,:)); 


%  adding  a  Gaussian  random  term  to  make  the  wind  speed  spread 
if  wind_veetor(l)~=0 

new_wind_vector(  1  )=wind_vector(  l)*normrnd(  1  ,abs(spread/max(wind_vector(  1 )))); 
else 

new_wind_vector(  1  )=0; 
end 
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if  wind_vector(2)~=0 

new_wind_vector(2)=wind_vector(2)*normrnd(l,abs(spread/max(wind_vector(2)))); 

else 

new_wind_vector(2)=0; 

end 

if  wind_vector(3)~=0 

new_wind_vector(3)=wind_vector(3)*normrnd(l,abs(spread/max(wind_vector(3)))); 

else 

new_wind_vector(3)=0; 

end 


for  i=l:number_size 

%  opening  of  data  file  registered  in  matrixpos. 
load(sprintf('size%d.maf,i)); 

%  total  number  of  drops 
drops=matrix_T_Round(  3  ,i ) ; 

%  transforming  to  rectangular  coordinates 
[x,y,z]=sph2cart(matrixpos(3,:),matrixpos(2,:),matrixpos(1 ,:)); 

diameter=matrix_T_Round(  1  ,i); 


% - 

%  Drop  fall  speed  depends  on  drop  diameter  in  the  z  axis, 

%  Atlas-Ulbrich  approximation  is  used 

%  Approximation  is  valid  in  the  diameter  range  5*e-4,  5*e-3 
% - 


wtmax=386.6*((matrix_T_Round(l,i)*0.001)A0.67); 

velocityz=(wtmax*(ones(l,drops)))*cos((rain_anglel*pi)/180); 

velocityx=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 

cos((rain_angle2*pi)/180); 

velocityy=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l  80); 

newest_wind_vector=[]; 


%  adding  Gaussian  random  term  to  make  the  wind  speed 
%  spread  between  different  drops  sizes 
if  new_wind_vector(l)~=0 

ne  west_wind_vector(  1 , :  )=new_wind_vector(  1 )  *normrnd(  1 , . . . 
abs(spread_2/max(new_wind_vector(l))),l, drops); 

else 

newest_wind_vector(l,:)=0; 

end 

if  new_wind_vector(2)~=0 

newest_wind_vector(2,:)=new_wind_vector(2)*normrnd(l,... 
abs(spread_2/max(new_wind_vector(2))),l, drops); 

else 

newest_wind_vector(2,:)=0; 
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end 

if  new_wind_vector(3)~=0 

newest_wind_vector(3,:)=new_wind_vector(3)*normrnd(l,... 
abs(spread_2/max(new_wind_vector(3))),l,  drops); 

else 

ne  west_wind_vector(  3 , :  )=0; 
end 

velocityxt=newest_wind_vector(  1 , :  )+velocityx; 
velocityyt=newest_wind_vector(2,:)+velocityy; 
velocityzt=newest_wind_vector(3,:)+velocityz; 

%  speed  registration 

m_velocity=[velocityxt;velocityyt;velocityzt]; 

save(sprintf('velocity%d.mat',i),'m_velocity'); 

deplacementx=velocityxt*(PRT_2-PRT_l); 

deplacementy=velocityyt*(PRT_2-PRT_l); 

deplacementz=velocityzt*(PRT_2-PRT_l); 

position_x=x+deplacementx; 

position_y=y+deplacementy; 

position_z=z+deplacementz; 

%  transformation  into  spherical  coordinates 

[theta,  phi,r]=cart2sph(position_x,position_y,position_z); 

%  computation  of  new  angles 
vector_azimuth=ones(l,drops)*(azimuth*pi)/180; 
vector_elevation=ones(l,drops)*(elevation*pi)/180; 
[xf,yf,zf]=sph2cart(vector_azimuth,vector_elevation,r); 

range=sqrt((position_x-xf).A2.+(position_y-yf).A2.+(position_z-zf).A2); 
%  approximation 
angle=atan(range ./ r) ; 

%  matrixpos 

matrixpos=[r;phi;theta;angle]; 

%  save  to  fde,  matrixpos 

s=sprintf(  'size%d.  mat' ,  i) ; 
save(s,  'matrixpos'); 


end 


%• 


%  add  to  E-field,  drop  and  rain  rate  comparison 
sum_E=[sum_E,total_E] ; 
weighted_E=[weighted_E,new_E]; 

end 


%• 
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% - 

%  Calculate  RCS  and  doppler 


%  calculate  RCS,  Z  and  Power  return  for  every  pulse  including  error 
running_Z=(Norm*abs(sum_E).A2)*wavelengthA4/piA5/abs(kw)A2/volume_resolution/... 
le-18+error_Z; 

power_return=Norm*abs(weighted_E).A2;  %  [W] 

%  calculate  average  Power  and  Z 
Z_avg=sum(running_Z)/p; 

Z_avg_plot=ones(  1  ,p)*Z_avg; 

power_return_avg=sum(Norm*abs(weighted_E).A2)/p;  %  [W] 

power_return_avg_plot=ones(l,p)*power_return_avg; 

%  calculate  doppler  frequency 
phase_vector_u=unwrap(phase_vector); 

doppler_fq=([phase_vector_u(2:2:2*p)-phase_vector_u(l:2:2*p-l)])/2/pi/PRT_l; 
number_vector=l  :p; 


% - 

%  Calculate  references 
%  Integral  form  of  Z 
ii=0: 100/1000: 100; 

zz=ii.A6.*8000.*exp(-4.1.*rain_rate.A(-.21).*ii); 

real_Z=ones(l,number_pulses)*10*logl0(trapz(ii,zz)); 

%  Adding  the  error  to  the  weighted  E  estimations 

error_power=error_rcs*P_t*GA2*wavelengthA2/piA2*(2*angle_3db*pi/180)A2*... 

(c*pulsewidth*le-6)/1024/log(2)/piA2/(distance*le3)A2/volume_resolution; 

%  Evaluate  average  Power  return  to  make  Z  estimation  including  error 
Z_estimation=power_return_avg/P_t/GA2/wavelengthA2*... 

piA2/(2*angle_3db*pi/180)A2/(c*pulsewidth*le-6)*1024*log(2)*(distance*le3)A2/... 

((kwA2)*(piA5)/(wavelengthA4))/le-18; 

Z_estimation_plot=ones(l,number_pulses)*10*logl0(Z_estimation+error_Z); 


% - 

%  Part  V,  Plot  of  results 
% - 


figure(l); 

subplot(2,l,l); 

plot(number_vector(  1  :p),doppler_fq(  1  :p),'-x',number_vector(  1  :p),ones(l,p)*... 
mean(dopplerfq),'— r'); 

title(['Doppler  frequency  for  all  pulse  pairs  for  PRF=  ',num2str(l/PRT_l),'  Hz']); 

xlabel('Pulse  pair'); 

ylabel('fd'); 

legend('Doppler  frequency  per  pulse  pair','Average  Doppler  estimation') 
grid  on 

subplot(2,l,2) 
hist(doppler_fq,  1 5) 
hold  on 

[ff,xx]  =  ksdensity(dopplerfq); 
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plot(xx,ff/max(ff)*max(hist(doppler_fq,15))) 
title('Doppler  Frequency  Histogram  with  Gaussian  Fit'); 
ylabel('Frequency'); 
xlabel('Doppler  Frequency  [Hz]'); 


figure  (2) 

subplot(2,l,l); 

%  also  converting  to  dBm 

plot(number_vector*PRT_2*  1 000, 1 0*log  1 0(power_retum+error_power)+3  0,'-xb',.. . 

number_vector*PRT_2*  1 000, 10*log  1 0(power_return_avg_plot+error_power)+30,'— r'); 
title(['Power  Return  for  all  ',num2str(p),'  pulses.  PRF=  ',num2str(l/PRT_2),'  Hz']); 
xlabel('Sample  time  [ms]'); 
ylabel('Power  [dBm]'); 

legend('Power  return',  'Average  Power  return', 0) 
grid  on 


subplot(2,l,2); 

plot(number_vector,  1 0*log  1 0(abs(running_Z)),'-xb',number_veetor,  10*... 
loglO(abs(Z_avg_plot)),'— r',number_vector,real_Z,'— k',... 
number_vector,Z_estimation_plot,'— b'); 
title('Estimated  Z  per  pulse'); 
xlabel('Sample  time'); 
ylabel('Z  [dBZ]'); 

legend('Simulation  Z,  (instantaneous)',  'Estimated  Z,  (average  return)',... 

'Z=\int  DA6N(D)dD', 'Estimated  Z,  Zeroth  moment', 0) 
grid  on 

figure(3) 

polarmatrix(  1:2:11 9)=zeros(  1,60); 
polarmatrix(2:2: 120)=phase_vector_u(  1 :60); 
lengthmatrix(  1:2:11 9)=zeros(  1 ,60); 
lengthmatrix(2 :2 : 1 2  0)=ones(  1,60); 
subplot(5,6,l) 

polar(polarmatrix(l:2),lengthmatrix(l:2),'-*') 
hold  on 

polar(polarmatrix(3 :4),lengthmatrix(3 :4),'— *r') 

title(['lA{st}.  v_r=  ',num2str(doppler_fq(l)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,2) 

polar(polarmatrix(5:6),lengthmatrix(5:6),'-*') 
hold  on 

polar(polarmatrix(7:8),lengthmatrix(7:8),'— *r') 

title(['2A{nd}.  v_r=  ',num2str(doppler_fq(2)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,3) 

polar(polarmatrix(9:10),lengthmatrix(9: 10),'-*') 
hold  on 

polar(polarmatrix(l  1:12),  lengthmatrix(  1 1:12),  ':*r') 

title(['3A{rd}.  v_r=  ',num2str(doppler_fq(3)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,4) 
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polar(polarmatrix(13:14),lengthmatrix(13:14),'-*') 
hold  on 

polar(polarmatrix(l  5 : 16),lengthmatrix(  1 5 : 1 6),':*r') 

title(['4A{thj .  v_r=  ',num2str(doppler_fq(4)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3 ) 

subplot(5,6,5) 

polar(polarmatrix(l  7 : 1 8),lengthmatrix(  1 7 : 18),'-*') 
hold  on 

polar(polarmatrix(l  9:20),  lengthmatrix(  1 9:20),  ':*r') 
titie(['5A{th}.  v_r=  ',num2str(doppler_fq(5)*wavelength/2),'  m/s']) 
hline  =  findobj(gca, 'Type', 'line'); 
set(hline,'LineWidth',3 ) 
subplot(5,6,6) 

polar(polarmatrix(2 1 :22),lengthmatrix(2 1 :22),'-*') 
hold  on 

polar(polarmatrix(23:24),lengthmatrix(23:24),':*r') 

title(['6A{thj .  v_r=  ',num2str(doppler_fq(6)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,7) 

polar(polarmatrix(25:26),lengthmatrix(25:26),'-*') 
hold  on 

polar(polarmatrix(27:28),lengthmatrix(27:28),':*r') 

title(['7A{th}.  v_r=  ',num2str(doppler_fq(7)*wavelength/2),'  m/s']) 

hline  =  findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,8) 

polar(polarmatrix(29:30),lengthmatrix(29:30),'-*') 
hold  on 

polar(polarmatrix(3 1 :32),lengthmatrix(3 1 :32),':*r') 

tide(['8A{th}.  v_r=  ',num2str(doppler_fq(8)*wavelength/2),'  m/s']) 

hline  =  findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,9) 

polar(polarmatrix(33:34),lengthmatrix(33:34),'-*') 
hold  on 

polar(polarmatrix(35:36),lengthmatrix(35:36),':*r') 

title(['9A{th| .  v_r=  ',num2str(doppler_fq(9)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,10) 

polar(polarmatrix(37:38),lengthmatrix(37:38),'-*') 
hold  on 

polar(polarmatrix(39:40),lengthmatrix(39:40),':*r') 

title(['10A{th}.  v_r=  ',num2str(doppler_fq(10)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,ll) 

polar(polarmatrix(4 1 :42),lengthmatrix(41 :42),'-*') 
hold  on 

polar(polarmatrix(43:44),lengthmatrix(43:44),':*r') 

title(['llA{th}.  v_r=  ',num2str(doppler_fq(ll)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,12) 
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polar(polarmatrix(45 :46),lengthmatrix(45 :46),'-*') 
hold  on 

polar(polarmatrix(47:48),lengthmatrix(47:48),':*r') 

title(['12A{th}.  v_r=  ’,num2str(doppler_fq(12)*wavelength/2),'  m/s']) 

hline  =  fmdobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3 ) 

subplot(5,6,13) 

polar(polarmatrix(49:50),lengthmatrix(49:50),'-*') 
hold  on 

polar(polarmatrix(51:52),lengthmatrix(51:52),':*r') 

title(['13A{th}.  v_r=  ',num2str(doppler_fq(13)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3 ) 

subplot(5,6,14) 

polar(polarmatrix(53:54),lengthmatrix(53:54),'-*') 
hold  on 

polar(polarmatrix(55:56),lengthmatrix(55:56),':*r') 

title(['14A{th}.  v_r=  ',num2str(doppler_fq(14)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,15) 

polar(polarmatrix(57:58),lengthmatrix(57:58),'-*') 
hold  on 

polar(polarmatrix(59:60),lengthmatrix(59:60),':*r') 

title(['15A{th}.  v_r=  ',num2str(doppler_fq(15)*wavelength/2),'  m/s']) 

hline  =  findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,16) 

polar(polarmatrix(61:62),lengthmatrix(61:62),'-*') 
hold  on 

polar(polarmatrix(63:64),lengthmatrix(63:64),':*r') 

title(['16A{th}.  v_r=  ',num2str(doppler_fq(16)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,17) 

polar(polarmatrix(65:66),lengthmatrix(65:66),'-*') 
hold  on 

polar(polarmatrix(67:68),lengthmatrix(67:68),':*r') 

title(['17A{th}.  v_r=  ',num2str(doppler_fq(17)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,18) 

polar(polarmatrix(69:70),lengthmatrix(69:70),'-*') 
hold  on 

polar(polarmatrix(71:72),lengthmatrix(71:72),':*r') 

title(['18A{th}.  v_r=  ',num2str(doppler_fq(18)*wavelength/2),'  m/s']) 

hline  =  findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,19) 

polar(polarmatrix(73:74),lengthmatrix(73:74),'-*') 
hold  on 

polar(polarmatrix(75:76),lengthmatrix(75:76),':*r') 

title(['19A{th}.  v_r=  ',num2str(doppler_fq(19)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,20) 
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polar(polarmatrix(77:78),lengthmatrix(77:78),'-*') 
hold  on 

polar(polarmatrix(79:80),lengthmatrix(79:80),':*r') 

title(['20A{th}.  v_r=  ',num2str(doppler_fq(20)*wavelength/2),'  m/s']) 

hline  =  fmdobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3 ) 

subplot(5,6,21) 

polar(polarmatrix(81:82),lengthmatrix(81:82),'-*') 
hold  on 

polar(polarmatrix(83:84),lengthmatrix(83:84),':*r') 

tide(['21A{th}.  v_r=  ',num2str(doppler_fq(21)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3 ) 

subplot(5,6,22) 

polar(polarmatrix(85:86),lengthmatrix(85:86),'-*') 
hold  on 

polar(polarmatrix(87:88),lengthmatrix(87:88),':*r') 

title(['22A{th}.  v_r=  ',num2str(doppler_fq(22)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,23) 

polar(polarmatrix(89:90),lengthmatrix(89:90),'-*') 
hold  on 

polar(polarmatrix(91:92),lengthmatrix(91:92),':*r') 

title(['23A{th}.  v_r=  ',num2str(doppler_fq(23)*wavelength/2),'  m/s']) 

hline  =  findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,24) 

polar(polarmatrix(93:94),lengthmatrix(93:94),'-*') 
hold  on 

polar(polarmatrix(95:96),lengthmatrix(95:96),':*r') 

title(['24A{th}.  v_r=  ',num2str(doppler_fq(24)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,25) 

polar(polarmatrix(97:98),lengthmatrix(97:98),'-*') 
hold  on 

polar(polarmatrix(99 : 1 00),lengthmatrix(99 : 1 00),' :  *r') 

title(['25A{th}.  v_r=  ',num2str(doppler_fq(25)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,26) 

polar(polarmatrix(  101:1 02),lengthmatrix(  101:1 02),'-*') 
hold  on 

polar(polarmatrix(  103:1 04),lengthmatrix(  103:1 04),':  *r') 

title(['26A{th}.  v_r=  ',num2str(doppler_fq(26)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,27) 

polar(polarmatrix(  105:1 06),lengthmatrix(  105:1 06),'-*') 
hold  on 

polar(polarmatrix(  107: 1 08),lengthmatrix(  1 07 : 1 08),':  *r') 

title(['27A{th}.  v_r=  ',num2str(doppler_fq(27)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,28) 
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polar(polarmatrix(  109:1 1 0),lengthmatrix(  1 09 : 1 1 0),'-*') 
hold  on 

polar(polarmatrix(  1 1 1:1 12),lengthmatrix(l  1 1 : 1 12),':*r') 

title(['28A{th}.  v_r=  ',num2str(doppler_fq(28)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,29) 

polar(polarmatrix(l  13:1 14),lengthmatrix(l  13:1 14),'-*') 
hold  on 

polar(polarmatrix(  115:1 1 6),lengthmatrix(  115:11 6),':  *r') 

title(['29A  {th} .  v_r=  ',num2str(doppler_fq(29)*wavelength/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,30) 

polar(polarmatrix(  117:11 8),lengthmatrix(  117:118),'-*') 
hold  on 

polar(polarmatrix(  119:1 20),lengthmatrix(  119:1 20),':  *r') 
title(['30A {th} .  v_r=  ',num2str(doppler_fq(30)*wavelength/2),'  m/s']) 
hline  =  findobj(gca,'Type','line'); 
set(hline,'LineWidth',3) 

mean(dopplerfq) 

std(dopplerfq) 

3)  WR_fq_step 

%  Weather  Radar  Signal  Simulator,  developed  version 
%  Frequency  step 
%  Ulf  Schroder 
%  August,  2005 

% - 


%  Main  program  for  The  NPS/Marielle  Gosset  Weather  Radar  Signal  Simulator. 


% - 

%  DESCRIPTION 
%  Part  I 

%  From  Radar  Data: 

%  -  Generate  Beam  Resolution  cell  approximation 

%  -  Generate  Cube  adding  a  margine  to  the  Beam  Resolution  cell  to  allow  for 
%  the  rain  drops  to  fall  trough  the  through  beam  over  time. 

% 

%  Part  II 

%  From  Rain  Parameters: 

%  Generating  the  fundamental  T  Matrix  concisting  of 
%  -  Diameter  Vector  and  Number  of  drops  in  each  diameter  interval 
%  using  the  Marshall-Palmer  exponential  approximation. 

%  -  the  total  number  of  drops  for  each  size  per  mA3 
%  The  backscatter  RCSs  for  each  Diameter  size  is  put  in  a  separate  RCS 
%  Matrix  to  account  for  the  differences  due  to  frequency 

%  Part  III 

%  Randomly  place  all  drops  in  cube 
%  Part  IV 


114 


%  For  every  pulse  the  coherent  Electric  field  return  is  calculated.  Between 
%  pulses  all  drops  are  moved. 

%  Part  V 

%  Calculate  and  plot  results 

% - 

clear  all 
close  all 
clc 

%  Input  Data 

c=3e8;  %  Speed  of  light  [m/s] 

%  Radar  Data 

P_t=  1 0e3 ;  %  T ransmitted  power 

G_dB=20;  %  Maximum  gain  in  dB 

G=10A(G_dB/10);  %  Maximum  gain 

angle_3db=0.20;  %  3  dB  Beam  angle  [degrees],  (half  half  power  beamwidth) 

pulsewidth=0.20;  %  Pulsewidth  [microsecond] 

elevation=3;  %  Elevation  angle  [degrees] 

margine=5;  %  margin  for  resolution  cell[m] 

azimuth=0;  %  Azimuth  angle  [rad] 

number_pulses=40;  %  Number  of  pulses 

PRT_l=l/5000;  %  Doppler  Pulse  Repetition  Time  [seconds] 

PRT_2=l/200;  %  Avg  Power  Pulse  Repetition  Time  [seconds] 

distance=0.5;  %  Range  from  the  radar  to  the  resolution  cell  [Km] 
wavelength=0.1;  %  Wavelength  [m] 

frequency=c/wavelength;  %  Frequency  [Hz] 

% - 

%Introduced  in  stepped  version 
fq_step=le6;  %  Frequency  step  in  [Hz] 

delay=50e-6;  %  Delay  for  changing  fq  in  [s] 

n_steps=3 ;  %  Number  of  Frequency  steps 

% - 

%  Rain  parameters 

diam_low=0.05;  %  Lowest  diameter  of  Precipitation  [mm] 

diam_lim=l;  %  Divider  between  low  and  mid  [mm] 

diam_mid=5;  %  Divider  between  mid  and  high  [mm] 

diam_high=15;  %  Highest  diameter  of  Precipitation  [mm] 
resolutionl=0.08;  %  Lower  diameter  resolution  [mm] 
resolution2=0.3;  %  Higher  diameter  resolution  [mm] 
resolution3=0.8;  %  Highest  diameter  resolution  [mm] 
sort_approx=l;  %  RCS  approximation  model;  [l]=Rayleigh 
kw=0.93;  %  abs(K)A2  approx  for  water 

rain_rate=50;  %  Rain  Rate  [mm/h] 
rain  angleM);  %  degrees  Rain  Angle  (Phi) 

rain_angle2=0;  %  degrees  Rain  Angle  (Theta) 

wind_vector=[-10;0;0];  %  Wind  speed  vector  (x,y,z)  [m/s] 
spread=0;  %  Wind  speed  spread  [m/s],  (all  directions) 

spread_2=4;  %  Drop  velocity  spread  [m/s], 

Norm- 1 00;  %  Normalization  coefficient 

%  Directional  cosines 
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u=sin(pi/2-elevation*pi/180)*cos(azimuth*pi/180); 

v=sin(pi/2-elevation*pi/180)*sin(azimuth*pi/180); 

w=cos(pi/2-elevation*pi/180); 


% - 

%  Part  I,  Generate  the  resolution  volume 

% - 

%  Approximation  of  the  resolution  cell  of  the  pulse  circular  box, 

%  for  any  azimuth  angle 

% - 

%  INPUT 
%  -  distance 
%  -  angle_3dB 
%  -  pulsewidth 
%  -  elevation 
%  -  margine 
%  -  azimuth 

% - 

%  OUTPUT 

%  -  volume  resolution:  computed  volume  of  the  resolution  cell 
%  -  volume_box:  volume  of  the  box  where  the  drops  are  put 
%  -  box:  the  8  rectangular  coordinates  of  the  box 
%  -  coord_vol_res:  the  coordinates  of  the  resolution  cell  where  the  drops 
%  are  taken. 

% - 


% - 

%  Creating  the  circular  beam  resolution  cell 
range_res=(c*pulsewidth*le-6)/2; 

%  Range  resolution  (ct/2)  [m] 

%  small  angle  approximation  allows  for  following 
angle_rad=(angle_3  db  *pi)/ 1 8  0 ; 

%  3  dB  beamwidth  [rad] 
radius  l=(tan(angle_rad)*(distance*  1000)); 

%  Radius  of  beam  at  R=distance 
radius2=(tan(angle_rad)*((distance*  1000)+range_res));  ... 

%  Radius  of  beam  at  R=distance  +  range  resolution 

surfacel=pi*radiuslA2; 

surface2=pi*radius2A2; 

%  volume  of  cone  =  base  surface  times  height  for  all  three 
vo  lume_cone  1  =(  surface  1  *(distance*  1 000))/3 ; 
volume_cone2=(surface2*((distance*1000)+range_res))/3; 
volume_resolution=(volume_cone2-volume_conel);  %  Beam  Cone  Volume 

%  computation  of  the  resolution  cell  coordinates,  spherical  coordinates 
%  elevation  in  radians 

elevation_rad=(elevation*pi)/180;  %  Elevation  [rad] 

azimuth_rad=(azimuth*pi)/180;  %  Azimuth  [rad] 

%  Matrix  used  when  comparing  wheather  a  drop  is  inside  or  outside 
%  beamwidth 

coord_vol_res=[distance*  1000, (distance*  1000)+range_res;elevation_rad-angle_rad,... 
elevation_rad+angle_rad;azimuth_rad-angle_rad,azimuth_rad+angle_rad]; 
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%• 


% - 

%  Creating  the  box 

%  computation  of  the  box  and  its  coordinates 

coin  l=[azimuth_rad+angle_rad,elevation_rad+angle_rad, distance*  1000]; 
coin2=[azimuth_rad-angle_rad,elevation_rad+angle_rad,  distance*  1000]; 
coin3=[azimuth_rad+angle_rad,elevation_rad-angle_rad,  distance*  1000]; 
coin4=[azimuth_rad-angle_rad,elevation_rad-angle_rad,distance*1000]; 
coin5=[azimiith_rad+angle_rad,elevation_rad+angle_rad,distance*1000+range_res]; 
coin6=[azimuth_rad-angle_rad, elevation_rad+angle_rad, distance*  1000+range_res]; 
coin7=[azimuth_rad+angle_rad,elevation_rad-angle_rad, distance*  1000+range_res]; 
coin8=[azimuth_rad-angle_rad,elevation_rad-angle_rad, distance*  1000+range_res]; 

[x  1  ,y  1  ,zl  ]=sph2cart(coin  1  ( 1  ),coin  1  (2), coin  1  (3 )); 
[x2,y2,z2]=sph2cart(coin2(l),coin2(2),coin2(3)); 

[x3  ,y3  ,z3]=sph2cart(coin3  ( 1  ),coin3  (2),coin3  (3 )); 
[x4,y4,z4]=sph2cart(coin4(l),coin4(2),coin4(3)); 
[x5,y5,z5]=sph2cart(coin5(l),coin5(2),coin5(3)); 
[x6,y6,z6]=sph2cart(coin6(l),coin6(2),coin6(3)); 
[x7,y7,z7]=sph2cart(coin7(l),coin7(2),coin7(3)); 
[x8,y8,z8]=sph2cart(coin8(l),coin8(2),coin8(3)); 

xmax=max([xl,x2,x3,x4,x5,x6,x7,x8]); 

xmin=min([xl,x2,x3,x4,x5,x6,x7,x8]); 

ymax=max([yl,y2,y3,y4,y5,y6,y7,y8]); 

ymin=min([y  1  ,y2,y3,y4,y5,y6,y7,y8]); 

zmax=max([zl,z2,z3,z4,z5,z6,z7,z8]); 

zmin=min([zl,z2,z3,z4,z5,z6,z7,z8]); 

%  the  box  with  its  margin 

coin_l=[xmin-margine,ymin-margine,zmin-margine]; 

coin_2=[xmin-margine,ymax+margine,zmin-margine]; 

coin_3=[xmax+margine,ymax+margine,zmin-margine]; 

coin_4=[xmax+margine,ymin-margine,zmin-margine]; 

coin_5=[xmin-margine,ymin-margine,zmax+margine]; 

coin_6=[xmin-margine,ymax+margine,zmax+margine]; 

coin_7=[xmax+margine,ymax+margine,zmax+margine]; 

coin_8=[xmax+margine,ymin-margine,zmax+margine]; 

box=[coin_r,coin_2',coin_3',coin_4',coin_5',coin_6',coin_7',coin_8']; 

%  calculate  the  volume  of  the  box 

volume_box=abs((xmax+margine-(xmin-margine))*(ymax+margine-(ymin-margine))... 

*(zmax+margine-(zmin-margine))); 


cube=[coin_r,coin_2',coin_3',coin_4',coin_r,coin_5',coin_6',coin_2',... 

coin_6',coin_7',coin_3',coin_7',coin_8',coin_4',coin_8',coin_5']; 

% - 

% - 

%  Part  II,  Generate  the  fundamental  T  Matrix 

% - 

%  This  part  produces  a  matrix  of  all  drops  that  contribute  to  the  total 
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%  reflectivity  of  the  resolution  cell 

% - 

%  First  all  sizes  are  created,  with  two  vectors  of  different  size  resolution. 

%  The  size  vectors  represents  all  diameters.  Next  the  backscatter  RCS  for  each 
%  size  is  computed. 

%  Using  the  Marshall-Palmer  exponential  approximation,  the  number  of  drops  of 
%  each  diameter  is  computed. 

% - 

%  INPUT 
%  -  diam  low 
%  -  diam_lim 
%  -  diam_high 
%  -  resolution  1 
%  -  resolution2 
%  -  sort_approx 
%  -  wavelength 
%  -  Kw 
%  -  rain_rate 

% - 

%  OUTPUT 
%  -  number_drops 
%  Matrix  T : 

%  -  first  row  all  sizes 
%  -  second  row  all  the  backscatter  RCSs 
%  -  third  row  the  total  number  of  drops  for  each  size  per  mA3 
%  the  columns  starting  from  the  smallest  size 

% - 

%  Create  dropsize  vector 

%  Precipitaion  Diameter  vector  (lowest<D<divider) 
vector_size_l=diam_low:resolutionl:diam_lim; 

%  Precipitaion  Diameter  vector  (divider<D<mid) 
vector_size_2=diam_lim:resolution2:diam_mid; 

%  Precipitaion  Diameter  vector  (mid<D<highest) 
vector_size_3=diam_mid:resolution3:diam_high; 

%  Precipitaion  Diameter  Vector 

vector_size=[vector_size_l,vector_size_2,vector_size_3]; 


% - 

%  Create  RCS  matrix  for  all  Rain  drops  of  Precipitation  Vector  and  for  all 
%  frequencies. 

fq_rr= frequency  :fq_step:frequency+(n_steps-l)*fq_step; 
rcs(l:n_steps,:)=((kwA2)*(piA5)/(cA4)*fq_rr'.A4)*((vector_size*le-3).A6); 
% - 


% - 

%  The  Marshall  Palmer  approximation  is  used  to  create  Precipitaion  Diameter 
%  Distribution  Vector  (Drop  Size  Distribution  N(D)  Vector) 

%  The  number  of  drops  are  scaled  down  by  a  factor  (Norm)  to  improve  the 
%  speed  of  the  simulator 

% - 

lambda=4.1*(rain_rateA(-0.21));  %[mm-l] 

number=8000*exp(-lambda*vector_size)/Norm; 
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vector_resl=ones(l,length(vector_size_l))*resolutionl; 

vector_res2=ones(l,length(vector_size_2))*resolution2; 

vector_res3=ones(l,length(vector_size_3))*resolution3; 

%  delta(Diameter)  dD 

vector_res=[vector_resl,vector_res2,vector_res3]; 

%  the  approximate  "true"  number  of  drops  N(D)dD  vector 
number=number.*vector_res; 

% - 


%■ 


%  create  matrix  T  without  res 
T=[vector_size;number]; 


%■ 


% - 

%  Change  matrix  T  to  matrix  'matrixTRound' 

% - 

%  This  part  computes  the  total  number  of  drops  by  utilizing  the  resolution 
%  cell  and  the  matrix_T_Round 

%  The  only  difference  between  the  matrix  T  Round  and  T,  is  the  fact  that 
%  matrix_T_Round  has  the  total  number  of  drops  on  the  third  row  and  they 
%  are  rounded  of  to  nearest  integer  value. 

% - 


matrix_T_Round=T ; 

matrix_T_Round(2,:)=round(matrix_T_Round(2,:)*volume_box); 

%  Keep  track  of  errors  due  to  round  off 
diff=T(2,:)*volume_box-matrix_T_Round(2,:); 
error_rcs=diff*rcs'*Norm/volume_box*volume_resolution; 
error_Z=error_rcs/volume_resolution*cA4./fq_rr.A4/piA5/abs(kw)A2/le-18; 


% - 

%  Part  III,  Initial  positioning  of  the  drops  within  the  volume 

% - 

%  This  part  generates  the  initial  position  of  the  drops  in  the  volume  of 
%  the  box. 

%  The  positions  are  recorded  in  the  file  "sizexx.maf  taking  into  account 
%  the  azimuth  angle 

% - 

%  INPUT 
%  -  box 

%  -  matrix_T_Round 
%  -  elevation 
%  -  azimuth 

% - 

%  OUTPUT 

%  All  files  including  the  sizes.  Each  contains  the  positions  of  all  drops 
%  for  a  given  diameter.  The  data  recoreded  as  follows: 

%  -  radial  distance  r 
%  -  elevation  phi 
%  -  azimuth  theta 

%  -  angle  with  respect  to  beam  center.  This  parameter  will  be 
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%  used  to  identify  drops  which  are  in  the  beam. 
% - 


%  Randomly  put  the  calculated  number  of  drops  of  each  diameter  in  the  box 
number_size=length(matrix_T_Round(l,:)); 

for  i=l:number_size 

%  select  number  of  drops  of  certain  Diameter 
drop  s=matrix_T_Round(2 ,  i ) ; 

%  create  as  many  positions  on  the  x  axis  as  #  drops  per  size 
elements_x=rand(  1  .drops); 

%  put  them  in  the  range  of  interest 
%  Random  placement  x  for  every  drop 

position_x=box(  1 , 1  )+(box(  1 ,3)-box(  1 , 1  ))*elements_x; 

%  create  as  many  positionson  the  y  axis  as  #  drops  per  size 
elements_y=rand(l,  drops); 

%  put  them  in  the  range  of  interest 
%  Random  placement  y  for  every  drop 

position_y=box(2, 1  )+(box(2,2)-box(2, 1  ))*elements_y ; 

%  create  as  many  positions  on  the  z  axis  as  #  drops  per  size 
elements_z=rand(  1  .drops); 

%  put  them  in  the  range  of  interest 
%  Random  placement  z  for  every  drop 

position_z=box(3,l)+(box(3,5)-box(3,l))*elements_z; 

%  matrix  of  results 

%  Creating  a  drop  Matrix  using  spherical  coordinates 

[theta,  phi,r]=cart2sph(position_x,position_y,position_z); 

%  Calculating  the  range,  phi,  theta  and  angle  to  every  drop  with 
%  reference  to  phi_0  and  theta_0.  Used  to  estimate  weighted  power  return. 
vector_azimuth=ones(l,drops)*((azimuth*pi)/180); 
vector_elevation=ones(  1  ,drops)*((elevation*pi)/ 180); 
[xf,yf,zfj=sph2cart(vector_azimuth,vector_elevation,r); 
range=sqrt((position_x-xf).A2.+(position_y-yf).A2.+(position_z-zf).A2); 

%  approximation 
angle=atan(range./r); 

%  Saving  the  drop  positions  in  size%d.mat 
%  matrixpos 

matrixpos=[r;phi;theta;angle]; 
save(sprintf('size%d. mat'd), 'matrixpos'); 
end 

% - 
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% - 

%  Part  IV,  Calculate  Coherent  Electric  field  return. 
% - 


%  For  every  pulse  the  coherent  Electric  field  return  is  calculatet.  Between 
%  pulses  all  drops  are  moved. 


% 


%  Run  the  simulation  for  specified  number  of  pulses 
% - 


%  Zero  all  vectors  and  matrices  that  are  going  to  be  used 
phase_vector=zeros(  1 ,2*number_pulses); 
imaginary_return=zeros(  1 ,2*number_pulses); 
sum_E=zeros(n_steps,ceil(number_pulses/n_steps)); 
weighted_E=zeros(n_steps,ceil(number_pulses/n_steps)); 
nmning_Z=[]; 

%  Initiate  reference  variables 

p=0;  %  Used  to  separate  pulses 

hh=0;  %  Used  to  separate  avg  power  samples 

while  p<number_pulses 

step_counter=0;  %  Used  to  separate  frequencies 

hh=hh+l; 

for  fq_new=frequency:fq_step:frequency+(n_steps-l)*fq_step 

k=2*pi*fq_new/c; 

p=p+l; 

step_counter=step_counter+ 1 ; 

disp(sprintf('calculate  the  power-return  of  pulse  %d\n',p)); 


% - 

%  computes  the  (weighted)  I  and  Q  return 

% - 

%  computes  the  total  number  of  drops 

% - 

%  INPUT 
%  -  Matrix  T 
%  -  wavelength 
%  -  angle_3dB 
%  -  distance 
%  -  coord_vol_res 
%  -  number_pulses 

% - 

%  OUTPUT 

%  -  total  E:  the  received  complex  Electric  field) 

%  -  power:  the  power 
%  -  voltage:  the  voltage 

%  -  drops_beam:  the  number  of  drops  in  the  beam 

% - 

%  Only  the  drops  which  are  in  the  resolution  cell  are  selected  for  the 
%  computation 

% - 
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number_size=length(matrix_T_Round(  1 

total_E=0; 

total_phase=0; 

new_E=0; 

for  i=l:number_size 

%  opening  the  file  of  registered  data  matrixpos. 
load(sprintf('size%d.mat',i)); 

%  removing  drops  that  are  outside  the  resolution  cell 

number_drops_per_size=length(matrixpos(l,:)); 

dropmatrix=[]; 

%  Filling  the  drop  matrix  with  all  drops  inside  the  beam 

%  resolution  volume 

for  d=l  :number_drops_per_size 

if  ((matrixpos(4,d)<=(angle_3db*pi/180))&(matrixpos(l,d)>=... 
coord_vol_res(  1 , 1  ))&(matrixpos(  1  ,d)<=coord_vol_res(  1,2))); 
dropmatrix=[dropmatrix,matrixpos( :  ,d)] ; 
end 
end 


% - 

%  Convert  to  cartesian  coordinates 
if  isempty(dropmatrix)==0 

[position_x,position_y,position_z]=sph2cart(dropmatrix(3 , :),... 
dropmatrix(2,  :),dropmatrix(  1, :)); 

%  Create  comparison  vectors 

vector_azimuth=ones(l,length(dropmatrix(l,:)))*((azimuth*pi)/180); 

vector_elevation=ones(l,length(dropmatrix(l,:)))*((elevation*pi)/180); 


%  Create  a  phase  vector  for  the  drops  of  current  Diameter 

r_phase=position_x.*u+position_y.*v+position_z.*w; 

phase=exp((2*j*k).*r_phase); 

%  Sum  the  contributions  to  the  E-field 

total_phase=total_phase+sum(phase); 

total_E=total_E+sum(sqrt(rcs(step_counter,i)).*phase); 


% - 

%  ceate  weighted  power  return,  concidering  range  (r),  theta,  phi 
%  with  reference  to  boresight  and  radar  parameters 
E_weighted=sqrt(P_t)*G*c/fq_rr(step_counter)/(sqrt(4*pi)A3)./... 
dropmatrix(l,:).A2.*(exp(-4*log(2).*((dropmatrix(3, :)-... 
vector_azimuth).A2/(2*angle_3db*pi/180)A2+((dropmatrix(2,:)... 
-vector_elevation).A2/(2*angle_3db*pi/180)A2)))); 
new_E=new_E+sum(E_weighted.*(sqrt(rcs(step_counter,i)).*phase)); 

end 


%• 
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phase_vector(2*p-l)=ANGLE(total_phase); 

imaginary_return(2*p-l)=total_phase; 

%  add  to  E-field  and  drop-  and  rate-vector 

sum_E(step_counter,hh)=total_E; 

weighted_E(step_counter,hh)=new_E; 


%■ 


% - 

%  This  part  simulates  the  movement  of  the  drops  and  gives  new  positions. 
%  First  using  the  Doppler  PRT  (PRT1) 

% - 

%  INPUT 

%  -  matrix_T_Round 
%  -  PRT 
%  -  azimuth 
%  -  elevation 
%  -  rain_anglel 
%  -  rain_angle2 
%  -  diam_lim 
%  -  resolution  1 
%  -  resolution2 

%  -  wind_vector 

% - 

%  OUTPUTS 

%  Output  files  containing  sizes  with  new  positions 
%  structured  the  same  way  as  in  previous 

%  there  is  also  an  output  file  for  speed  (one  for  each  size)  for  the  pulse 
%  pair. 

% - 


number_size=length(matrix_T_Round(l,:)); 

%  adding  a  Gaussian  random  term  to  make  the  wind  speed  spread 
if  wind_vector(l)~=0 

new_wind_vector(  1  )=wind_vector(  1  )*normrnd(  l,abs(spread/max(wind_vector(  1 )))); 
else 

new_wind_vector(  1  )=0; 
end 

if  wind_vector(2)~=0 

new_wind_vector(2)=wind_vector(2)*normrnd(l,abs(spread/max(wind_vector(2)))); 

else 

new_wind_vector(2)=0; 

end 

if  wind_vector(3)~=0 

new_wind_vector(3)=wind_vector(3)*normrnd(l,abs(spread/max(wind_vector(3)))); 

else 

new_wind_vector(3  )=0; 
end 

for  i=l:number_size 

%  opening  of  data  file  registered  in  matrixpos. 
load(sprintf('size%d.maf,i)); 
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%  total  number  of  drops 
drop  s=matrix_T_Round(2 ,  i ) ; 

%  transforming  to  rectangular  coordinates 
[x,y,z]=sph2cart(matrixpos(3,:),matrixpos(2,:),matrixpos(1 

diameter=matrix_T_Round(  1  ,i); 


% - 

%  Drop  fall  speed  depends  on  drop  diameter  in  the  z  axis, 

%  Atlas-Ulbrich  approximation  is  used 

%  Approximation  is  valid  in  the  diameter  range  5*e-4,  5*e-3 
% - 


wtmax=386.6*((matrix_T_Round(l,i)*0.001)A0.67); 

velocity  z=(wtmax*(ones(l,drops)))*cos((rain_anglel*pi)/l  80); 
velocityx=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
cos((rain_angle2*pi)/180); 

velocityy=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l  80); 


% - 

% - 

%  %  adding  another  Gaussian  random  term  to  make  the  wind  speed 

%  %  spread  between  different  drops  sizes 

%  if  new_wind_vector(  1  )~=0 

%  newest_wind_vector(  1  )=new_wind_vector(  1  )*normrnd(  1  ,abs(spread_2/... 

%  max(wind_vector(l)))); 

%  else 

%  newest_wind_veetor(l)=0; 

%  end 

%  if  new_wind_vector(2)~=0 

%  newest_wind_vector(2)=new_wind_vector(2)*normrnd(l,abs(spread_2/... 

%  max(wind_vector(2)))); 

%  else 

%  newest_wind_vector(2)=0; 

%  end 

%  if  new_wind_vector(3)~=0 

%  newest_wind_vector(  3  )=new_wind_vector(  3 )  *normrnd(  1  ,abs(  spread_2/. . . 

%  max(wind_vector(3)))); 

%  else 

%  newest_wind_vector(3)=0; 

%  end 

% 

%  %  If  spread  within  pulse 

%  velocityxt=newest_wind_vector(  l)-velocityx; 

%  velocityyt=newest_wind_vector(2)-velocityy; 

%  velocityzt=newest_wind_vector(3)-velocityz; 

% 

%  %  speed  registration 

%  m_velocity=[velocityxt;velocityyt;velocityzt]; 

%  save(sprintf('velocity%d.maf,i),'m_velocity'); 


%• 
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% - 

%  If  spread  for  all  drops  regardless  of  size 
newest_wind_vector=[]; 


%  adding  Gaussian  random  term  to  make  the  wind  speed 
%  spread  between  drops 
if  new_wind_vector(  1  )~=0 

newest_wind_vector(  1 , :  )=new_wind_vector(  1 )  * . . . 

normrnd(  1  ,abs(spread_2/max(new_wind_vector(  1 ))),  1 ,  drops); 

else 

newest_wind_vector(l,:)=0; 

end 

if  new_wind_vector(2)~=0 

newest_wind_vector(2,  :)=new_wind_vector(2)* . . . 

normrnd(  1  ,abs(spread_2/max(new_new_wind_vector(2))),  1 ,  drops); 

else 

newest_wind_vector(2,:)=0; 

end 

if  new_wind_vector(3)~=0 
newest_wind_vector(  3 , :  )=new_wind_vector(  3 )  * . . . 

normrnd(  1  ,abs(spread_2/max(new_wind_vector(3))),  1 ,  drops); 

else 

newest_wind_vector(  3 , :  )=0; 
end 

%  If  spread  within  pulse 
velocityxt=newest_wind_vector(l,:)-velocityx; 
velocityyt=newest_wind_vector(2,:)-velocityy; 
velocityzt=newest_wind_vector(3 , :  )-velocityz; 


%■ 

%■ 


deplacementx=velocityxt*PRT_l ; 
deplacementy=velocityyt*PRT_l ; 
deplacementz=velocityzt*PRT_l ; 

position_x=x+deplacementx; 

position_y=y+deplacementy; 

position_z=z+deplacementz; 

%  transformation  into  spherical  coordinates 

[theta,  phi,r]=cart2sph(position_x,position_y,position_z); 

%  computation  of  new  angles 
vector_azimuth=ones(l,drops)*(azimuth*pi)/180; 
vector_elevation=ones(l,drops)*(elevation*pi)/180; 
[xf,yf,zf]=sph2cart(vector_azimuth,vector_elevation,r); 

range=sqrt((position_x-xf).A2.+(position_y-yf).A2.+(position_z-zf).A2); 

%  approximation 

angle=atan(range./r); 

%  matrixpos 

matrixpos=[r;phi;theta;angle]; 
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%  save  to  file,  matrixpos 

s=sprintf('size%d.mat',i); 
save(s,  'matrixpos'); 


end 


% - 

%  Calculate  phase  for  pulse  pair 
total_phase=0; 

for  i=l:number_size 

%  opening  the  file  of  registered  data  matrixpos. 
load(sprintf('size%d.mat',i)); 

%  removing  drops  that  are  outside  the  resolution  cell 
number_drops_per_size=length(matrixpos(l,:)); 
dropmatrix=[  |; 

%  Filling  the  drop  matrix  with  all  drops  inside  the  beam 

%  resolution  volume 

for  d=l  :number_drops_per_size 

if  ((matrixpos(4,d)<=(angle_3db*pi/180))&(matrixpos(l,d)>=... 
coord_vol_res(  1 , 1  ))&(matrixpos(  1  ,d)<=coord_vol_res(  1,2))); 
dropmatrix=[dropmatrix,matrixpos( :  ,d)] ; 
end 
end 


% - 

%  Convert  to  cartesian  coordinates 
if  isempty(dropmatrix)==0 

[position_x,position_y,position_z]=sph2cart(dropmatrix(3 , :),... 
dropmatrix(2,  :),dropmatrix(  1, :)); 

%  Create  comparison  vectors 

vector_azimuth=ones(l,length(dropmatrix(l,:)))*((azimuth*pi)/180); 

vector_elevation=ones(l,length(dropmatrix(l,:)))*((elevation*pi)/180); 


%  Create  a  phase  vector  for  the  drops  of  current  Diameter 

r_phase=position_x.*u+position_y.*v+position_z.*w; 

phase=exp((2*j*k).*r_phase); 

%  Sum  the  contributions  to  the  E-field 
total_phase=total_phase+sum(phase); 
end 
end 

phase_vector(2*p)=ANGLE(total_phase); 
imaginary  _return(2*p)=total_phase; 


% - 

%  This  part  simulates  the  movement  of  the  drops  and  gives  new  positions. 
%  Using  the  doppler  PRT  and  the  delay  (PRT  2,  delay) 
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%• 


%■ 


number_size=length(matrix_T_Round(l,:)); 

%  adding  a  Gaussian  random  term  to  make  the  wind  speed  spread 
if  wind_veetor(l)~=0 

new_wind_vector(  1  )=wind_vector(  1  )*normrnd(  1  ,abs(spread/. . . 
max(wind_vector(l)))); 

else 

new_wind_vector(  1  )=0; 
end 

if  wind_vector(2)~=0 

new_wind_vector(2)=wind_vector(2)*normrnd(l,abs(spread/... 

max(wind_vector(2)))); 

else 

new_wind_vector(2)=0; 

end 

if  wind_vector(3)~=0 

new_wind_vector(3  )=wind_vector(3  )*normrnd(  1  ,abs(spread/. . . 
max(wind_vector(3)))); 

else 

new_wind_vector(3)=0; 

end 

for  i=l  :number_size 

%  opening  of  data  file  registered  in  matrixpos. 
load(sprintf('size%d. mat'd)); 

%  total  number  of  drops 
drop  s=matrix_T_Round(2  ,i ) ; 

%  transforming  to  rectangular  coordinates 
[x,y,z]=sph2cart(matrixpos(3 , :  ),matrixpos(2, :  ),matrixpos(  1 ,:)); 

diameter=matrix_T_Round(  1  ,i); 


% - 

%  Drop  fall  speed  depends  on  drop  diameter  in  the  z  axis, 

%  Atlas-Ulbrich  approximation  is  used 

%  Approximation  is  valid  in  the  diameter  range  5*e-4,  5*e-3 
% - 


wtmax=3  8  6 . 6  *  ((matrix_T_Round(  1  d)  *  0 . 00 1 )  A0. 6  7) ; 

velocityz=(wtmax*(ones(l,drops)))*cos((rain_anglel*pi)/180); 

velocityx=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 

cos((rain_angle2*pi)/180); 

velocityy=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l  80); 


%■ 

Vo- 
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%  %  adding  another  Gaussian  random  term  to  make  the  wind  speed 

%  %  spread  between  different  drops  sizes 

%  if  new_wind_vector(  1  )~=0 

%  newest_wind_vector(  1  )=new_wind_vector(  1  )*normrnd(  1  ,abs(spread_2/... 

%  max(wind_vector(l)))); 

%  else 

%  newest_wind_vector(l)=0; 

%  end 

%  if  new_wind_vector(2)~=0 

%  newest_wind_vector(2)=new_wind_vector(2)*normrnd(l,abs(spread_2/... 

%  max(wind_vector(2)))); 

%  else 

%  newest_wind_vector(2)=0; 

%  end 

%  if  new_wind_vector(3)~=0 

%  newest_wind_vector(3)=new_wind_vector(3)*normrnd(l,abs(spread_2/... 

%  max(wind_vector(3)))); 

%  else 

%  newest_wind_vector(3)=0; 

%  end 

% 

%  %  If  spread  within  pulse 

%  velocityxt=newest_wind_vector(  1  )-velocityx; 

%  velocityyt=newest_wind_vector(2)-velocityy; 

%  velocityzt=newest_wind_vector(3)-velocityz; 

% 

% - 

% - 

%  If  spread  for  all  drops  regardless  of  size 
newest_wind_vector=[] ; 

%  adding  Gaussian  random  term  to  make  the  wind  speed 
%  spread  between  drops 
if  new_wind_vector(  1  )~=0 

newest_wind_vector(  1 , :  )=new_wind_vector(  1 )  * . . . 

normrnd(  1  ,abs(spread_2/max(new_wind_vector(  1  ))),1 ,  drops); 

else 

newest_wind_vector(l,:)=0; 

end 

if  new_wind_vector(2)~=0 

newest_wind_vector(2,  :)=new_wind_vector(2)* . . . 

normrnd(  1  ,abs(spread_2/max(new_new_wind_vector(2))),  1 ,  drops); 

else 

newest_wind_vector(2,:)=0; 

end 

if  new_wind_vector(3)~=0 

newest_wind_vector(  3 , :  )=new_wind_vector(  3 )  * . . . 
normrnd(l,abs(spread_2/max(new_wind_vector(3))),l,  drops); 

else 

newest_wind_vector(  3 , :  )=0; 
end 
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%  If  spread  within  pulse 
velocityxt=newest_wind_vector(l,:)-velocityx; 
velocityyt=newest_wind_vector(2,:)-velocityy; 
velocityzt=newest_wind_vector(3,:)-velocityz; 


o/o- 

o/o- 


%  speed  registration 

m_velocity=[  velocityxt;velocityyt;velocityzt]; 
save(sprintf('velocity%d.mat',i),'m_velocity'); 

deplacementx=velocityxt*(PRT_l+delay); 

deplacementy=velocityyt*(PRT_l+delay); 

deplacementz=velocityzt*(PRT_l+delay); 

position_x=x+deplacementx; 

position_y=y+deplacementy; 

position_z=z+deplacementz; 

%  transformation  into  spherical  coordinates 

[theta,  phi,r]=cart2sph(position_x,position_y,position_z); 

%  computation  of  new  angles 
vector_azimuth=ones(l,drops)*(azimuth*pi)/180; 
vector_elevation=ones(l,drops)*(elevation*pi)/180; 
[xf,yf,zf]=sph2cart(vector_azimuth,vector_elevation,r); 

range=sqrt((position_x-xf).A2.+(position_y-yf).A2.+(position_z-zf).A2); 

%  approximation 

angle=atan(range./r); 

%  matrixpos 

matrixpos=[r;phi;theta;angle]; 

%  save  to  fde,  matrixpos 

s=sprintf(  'size%d.  mat' ,  i) ; 
save(s,  'matrixpos'); 


end 


%• 


% - 

%  This  part  simulates  the  movement  of  the  drops  and  gives  new  positions. 
%  Using  the  Avg  Power  PRT  reduced  by  the  sample  times  for  the  previous 
%  samples 

% - 


% - 

number_size=length(matrix_T_Round(  1 

%  adding  a  Gaussian  random  term  to  make  the  wind  speed  spread 
if  wind_vector(l)~=0 


129 


new_wind_vector(l)=wind_vector(l)*normrnd(l,abs(spread/... 

max(wind_vector(l)))); 

else 

new_wind_vector(  1  )=0; 
end 

if  wind_vector(2)~=0 

new_wind_vector(2)=wind_vector(2)*normrnd(l,abs(spread/... 

max(wind_vector(2)))); 

else 

new_wind_vector(2)=0; 

end 

if  wind_vector(3)~=0 

new_wind_vector(3)=wind_vector(3)*normrnd(l,abs(spread/... 

max(wind_vector(3)))); 

else 

new_wind_vector(3  )=0; 
end 


for  i=l:number_size 

%  opening  of  data  file  registered  in  matrixpos. 
load(sprintf('size%d.mat',i)); 

%  total  number  of  drops 
drop  s=matrix_T_Round(  2 ,  i) ; 

%  transforming  to  rectangular  coordinates 
[x,y,z]=sph2cart(matrixpos(3,:),matrixpos(2,:),matrixpos(l 

diameter=matrix_T_Round(  1  ,i); 


% - 

%  Drop  fall  speed  depends  on  drop  diameter  in  the  z  axis, 

%  Atlas-Ulbrich  approximation  is  used 

%  Approximation  is  valid  in  the  diameter  range  5*e-4,  5*e-3 
% - 


wtmax=386.6*((matrix_T_Round(l,i)*0.001)A0.67); 

velocity  z=(wtmax*(ones(l,drops)))*cos((rain_anglel*pi)/l  80); 
velocityx=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
cos((rain_angle2*pi)/180); 

velocityy=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l  80); 


% - 

% - 

%  %  adding  another  Gaussian  random  term  to  make  the  wind  speed 

%  %  spread  between  different  drops  sizes 

%  if  new_wind_vector(  1  )~=0 

%  newest_wind_vector(  1  )=new_wind_vector(  1  )*normrnd(  1  ,abs(spread_2/... 

%  max(wind_vector(l)))); 

%  else 

%  newest_wind_vector(l)=0; 

%  end 
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%  if  new_wind_vector(2)~=0 

%  newest_wind_vector(2)=new_wind_vector(2)*normrnd(l,abs(spread_2/. 

%  max(wind_vector(2)))); 

%  else 

%  newest_wind_vector(2)=0; 

%  end 

%  if  new_wind_vector(3)~=0 

%  newest_wind_vector(3)=new_wind_vector(3)*normrnd(l,abs(spread_2/. 

%  max(wind_vector(3)))); 

%  else 

%  newest_wind_vector(3)=0; 

%  end 

% 

%  %  If  spread  within  pulse 

%  velocityxt=newest_wind_vector(  l)-velocityx; 

%  velocityyt=newest_wind_vector(2)-velocityy; 

%  velocityzt=newest_wind_vector(3)-velocityz; 


% - 

% - 

%  If  spread  for  all  drops  regardless  of  size 
newest_wind_vector=[]; 


%  adding  Gaussian  random  term  to  make  the  wind  speed 
%  spread  between  drops 
if  new_wind_vector(l)~=0 

newest_wind_vector(  1 , :  )=ne  w_wind_vector(  1 )  * . . . 

normrnd(  1  ,abs(spread_2/max(new_wind_vector(  1 ))),  1 ,  drops); 

else 

newest_wind_vector(l,:)=0; 

end 

if  new_wind_vector(2)~=0 

newest_wind_vector(2,  :)=new_wind_vector(2)* . . . 
normrnd(l,abs(spread_2/max(new_new_wind_vector(2))),l,  drops); 

else 

newest_wind_vector(2,:)=0; 

end 

if  new_wind_vector(3)~=0 
newest_wind_vector(  3 , :  )=new_wind_vector(  3 )  * . . . 

normrnd(  1  ,abs(spread_2/max(new_wind_vector(3))),  1 ,  drops); 

else 

newest_wind_vector(  3 , :  )=0; 
end 

%  If  spread  within  pulse 
velocityxt=newest_wind_vector(l,:)-velocityx; 
velocityyt=newest_wind_vector(2,:)-velocityy; 
velocityzt=newest_wind_vector(3 , :  )-velocityz; 


%• 

%• 


%  speed  registration 

m_velocity=[velocityxt;velocityyt;velocityzt]; 
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save(sprintf('velocity%d.mat',i),'m_velocity'); 

deplacementx=velocityxt*(PRT_2-n_steps*(PRT_l+delay)); 

deplacementy=velocityyt*(PRT_2-n_steps*(PRT_l+delay)); 

deplacementz=velocityzt*(PRT_2-n_steps*(PRT_l+delay)); 

position_x=x+deplacementx; 

position_y=y+deplacementy; 

position_z=z+deplacementz; 

%  transfonnation  into  spherical  coordinates 
[theta,phi,r]=cart2sph(position_x,position_y,position_z); 

%  computation  of  new  angles 
vector_azimuth=ones(l,drops)*(azimuth*pi)/180; 
vector_elevation=ones(l,drops)*(elevation*pi)/180; 
[xf,yf,zf]=sph2cart(vector_azimuth,vector_elevation,r); 

range=sqrt((position_x-xf).A2.+(position_y-yf).A2.+(position_z-zf).A2); 

%  approximation 

angle=atan(range./r); 

%  matrixpos 

matrixpos=[r;phi;theta;angle]; 

%  save  to  fde,  matrixpos 

s=sprintf('size%d. mat',i); 
save(s,  'matrixpos'); 
end 

end 

end 


%■ 


% - 

%  Calculate  RCS  and  doppler 
% - 


% - 

%  calculate  RCS,  Z  and  Power  return  for  every  pulse  including  error 

% - 

for  o=l:n_steps 

running_Z(o,:)=(Norm*abs(sum_E(o,:))  A2)*cA4/fq_rr(o)A4/piA5/... 
abs(kw)A2/volume_resolution/ 1  e- 1 8+error_Z(o); 
end 

mnning_Z=reshape(running_Z,l,[]);  %  Transform  into  a  vector 

% - 

%  Adding  the  error  to  the  weighted  E  estimations 
% - 


error_power=error_rcs*P_t*GA2*cA2./fq_rr.A2/piA2*(2*angle_3db*pi/180)A2*... 

(c*pulsewidth*le-6)... 

/1 024/log(2)/piA2/(distance*  1  e3)A2; 
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for  o=l:n_steps 

power_return_matrix(o,:)=Norm*abs(weighted_E(o,:)).A2+error_power(o); 

end 

power_return=reshape(power_return_matrix, !,[]);%  Transform  into  a  vector 


% - 

%  calculate  average  Power  and  Z 

% - 

Z_avg=sum(running_Z)/p; 

Z_avg_plot=ones(  1  ,p)*Z_avg; 

power_return_avg=sum(sum(Norm*abs(weighted_E).A2,2)+error_power')/p;  %  [W] 
power_return_avg_plot=ones(l,p)*power_return_avg; 


% - 

%  calculate  doppler  frequency 
%  use  unwrap  to  make  sure  no  alias  is  presented 

% - 

phase_vector_u=unwrap(phase_vector); 

doppler_fq=([phase_vector_u(2:2:2*p)-phase_vector_u(l:2:2*p-l)])/2/pi/PRT_l; 
number_vector=l  :p; 


% - 

%  Calculate  references 
%  Integral  form  of  Z 

% - 

ii=0: 100/1000: 100; 

zz=ii.A6.*8000.*exp(-4.1.*rain_rate.A(-.21).*ii); 

real_Z=ones(l,p)*10*logl0(trapz(ii,zz)); 


% - 

%  Evaluate  average  Power  return  to  make  Z  estimation  including  error 

% - 

%  To  make  the  estimate  with  the  correct  frequency,  the  power  return  Matrix 
%  is  used.  For  the  plot  the  result  must  me  summed  and  divided  by  the  number 
%  of  samples  used  to  average 

Z_estimation=power_retum_matrix'*(fq_rr.A2)'/P_t/GA2/cA2*„. 

piA2/(2*angle_3db*pi/180)A2/(c*pulsewidth*le-6)*1024*log(2)*(distance*le3)A2/... 

((kwA2)*(piA5)/(wavelengthA4))/le-18; 

Z_estimation_plot=ones(l,p)*10*logl0((sum(Z_estimation)+sum(error_Z))/p); 


% - 

%  Part  V,  Plot  of  results 
% - 


figure(l); 

subplot(2,l,l); 

plot(number_vector(  1  :p),doppler_fq(  1  :p),'-x',number_vector(  1  :p),ones(l,p)*... 
mean(dopplerfq),'— r'); 

title(['Doppler  frequency  for  all  pulse  pairs  for  PRF=  ',num2str(l/PRT_l),'  Hz']); 

xlabel('Pulse  pair'); 

ylabel('fd'); 
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legend('Doppler  frequency  per  pulse  pairVAverage  Doppler  estimation') 
grid  on 


subplot(2,l,2) 

%  To  account  for  the  differences  in  wavelengths  a  correct  wavelength  vector 
%  is  created  and  multiplied  with  the  doppler  frequency  vector 
u=ones(  1  ,p/n_steps); 

lambida=[u*c/fq_rr(l);u*c/fq_rr(2);u*c/fq_rr(3)]; 
wave_l=reshape(lambida,l,[]); 
v_r=doppler_fq.  *  wave_l/2 ; 

hist(v_r,15) 
hold  on 

[ff,xx]  =  ksdensity(vr); 

plot(xx,ff/max(ff)*max(hist(v_r,15))) 

title('Radial  Velocity  Histogram  with  Gaussian  Fit'); 

xlabel('Radial  velocity  [m/s]'); 

ylabel('Frequency'); 

grid  on 


figure(2) 

subplot(2,l,l); 

kl=0:p/n_steps-l; 

power_number_vector=[kl*PRT_2;  kl*PRT_2+PRT_l+PRT_l+delay;  kl*PRT_2+2*... 
(PRT_  1  +PRT_  1  +delay )] ; 

power_number_vector=reshape(power_number_vector,l,[]); 

text_string=['Power  Return  for  all  ',num2str(p),'  pulses.','  Frequency  stepped  from 
num2str(fq_rr(l)/le9),'  GHz  to  ',num2str(fq_rr(n_steps)/le9),'  GHz  with 
num2str(fq_step/le6),'  MHz  steps.']; 

%  also  converting  to  dBm 

plot(power_number_vector*1000,10*logl0(power_return)+30,'-xb',... 

power_number_vector*1000,10*logl0(power_return_avg_plot)+30,'— r'); 
title(text_string); 
xlabel('Sample  time  [ms]'); 
ylabel('Power  [dBm]'); 

legend('Power  return',  'Average  Power  return', 0) 
grid  on 

subplot(2,l,2); 

plot(power_number_vector*1000,10*logl0(abs(iunning_Z)),'-xb',power_number_vector*... 
1000, 10*logl0(abs(Z_avg_plot)), '-.r',power_number_vector*  1000, real_Z,'-k',... 
power_number_vector*1000,Z_estimation_plot,'— b'); 
title('Estimated  Z  per  pulse'); 
xlabel('Sample  time'); 
ylabel('Z  [dBZ]'); 

legend('Simulation  Z,  (instantaneous)',  'Estimated  Z,  (average  return)',... 

'Z=\int  DA6N(D)dD', 'Estimated  Z,  Zeroth  moment',0) 
grid  on 

figure  (3) 

%  This  plot  only  works  if  the  number  of  pulse  pairs  exceed  30 
polarmatrix(  1:2:1 19)=zeros(l,60); 
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polarmatrix(2:2: 120)=phase_vector_u(  1 :60); 
lengthmatrix(  1:2:11 9)=zeros(  1 ,60); 
lengthmatrix(2 :2 : 1 2  0)=ones(  1,60); 
subplot(5,6,l) 

polar(polarmatrix(l:2),lengthmatrix(l:2),'-*') 
hold  on 

polar(polarmatrix(3 :4),lengthmatrix(3 :4),' — *r') 

title(['lA{st}.  v_r=  ',num2str(doppler_fq(l)*c/fq_rr(l)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3 ) 

subplot(5,6,2) 

polar(polarmatrix(5:6),lengthmatrix(5:6),'-*') 
hold  on 

polar(polarmatrix(7:8),lengthmatrix(7:8),'— *r') 

title(['2A{nd}.  v_r=  ',num2str(doppler_fq(2)*c/fq_rr(2)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,3) 

polar(polarmatrix(  9 : 1 0 ),  lengthmatrix(  9:10),'-*') 
hold  on 

polar(polarmatrix(l  1:12),  lengthmatrix(  1 1:12),  ':*r') 

title(['3A{rd}.  v_r=  ',num2str(doppler_fq(3)*c/fq_rr(3)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,4) 

polar(polarmatrix(13:14),lengthmatrix(13:14),'-*') 
hold  on 

polar(polarmatrix(l  5 : 16),lengthmatrix(  1 5: 1 6),':*r') 

title(['4A{th] .  v_r=  ',num2str(doppler_fq(4)*c/fq_rr(l)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,5) 

polar(polarmatrix(l  7 : 1 8),lengthmatrix(  1 7 : 18),'-*') 
hold  on 

polar(polarmatrix(l  9:20),  lengthmatrix(  1 9:20),  ':*r') 

title(['5A{th] .  v_r=  ',num2str(doppler_fq(5)*c/fq_rr(2)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,6) 

polar(polarmatrix(21 :22),lengthmatrix(21 :22),'-*') 
hold  on 

polar(polarmatrix(23:24),lengthmatrix(23:24),':*r') 

title(['6A{th| .  v_r=  ',num2str(doppler_fq(6)*c/fq_rr(3)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,7) 

polar(polarmatrix(25:26),lengthmatrix(25:26),'-*') 
hold  on 

polar(polarmatrix(27:28),lengthmatrix(27:28),':*r') 

title(['7A{thj .  v_r=  ',num2str(doppler_fq(7)*c/fq_rr(l)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,8) 

polar(polarmatrix(29:30),lengthmatrix(29:30),'-*') 
hold  on 

polar(polarmatrix(3 1 :32),lengthmatrix(3 1 :32),':*r') 
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title(['8A{th}.  v_r=  ',num2str(doppler_fq(8)*c/fq_rr(2)/2),'  m/s']) 
hline  =  findobj(gca, 'Type', 'line'); 
set(hline,'LineWidth',3 ) 
subplot(5,6,9) 

polar(polarmatrix(33:34),lengthmatrix(33:34),'-*') 
hold  on 

polar(polarmatrix(35:36),lengthmatrix(35:36),':*r') 
title(['9A{th}.  v_r=  ',num2str(doppler_fq(9)*c/fq_rr(3)/2),'  m/s']) 
hline  =  findobj(gca, 'Type', 'line'); 
set(hline,'LineWidth',3 ) 
subplot(5,6,10) 

polar(polarmatrix(37:38),lengthmatrix(37:38),'-*') 
hold  on 

polar(polarmatrix(39:40),lengthmatrix(39:40),':*r') 

title(['10A{th}.  v_r=  ',num2str(doppler_fq(10)*c/fq_rr(l)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,ll) 

polar(polarmatrix(41 :42),lengthmatrix(41 :42),'-*') 
hold  on 

polar(polarmatrix(43:44),lengthmatrix(43:44),':*r') 

title(['llA{th}.  v_r=  ',num2str(doppler_fq(ll)*c/fq_rr(2)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,12) 

polar(polarmatrix(45 :46),lengthmatrix(45 :46),'-*') 
hold  on 

polar(polarmatrix(47:48),lengthmatrix(47:48),':*r') 

title(['12A{th}.  v_r=  ',num2str(doppler_fq(12)*c/fq_rr(3)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,13) 

polar(polarmatrix(49:50),lengthmatrix(49:50),'-*') 
hold  on 

polar(polarmatrix(51:52),lengthmatrix(51:52),':*r') 

title(['13A{th}.  v_r=  ',num2str(doppler_fq(13)*c/fq_rr(l)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,14) 

polar(polarmatrix(53:54),lengthmatrix(53:54),'-*') 
hold  on 

polar(polarmatrix(55:56),lengthmatrix(55:56),':*r') 

title(['14A{th}.  v_r=  ',num2str(doppler_fq(14)*c/fq_rr(2)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,15) 

polar(polarmatrix(57:58),lengthmatrix(57:58),'-*') 
hold  on 

polar(polarmatrix(59:60),lengthmatrix(59:60),':*r') 

title(['15A{th}.  v_r=  ',num2str(doppler_fq(15)*c/fq_rr(3)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,16) 

polar(polarmatrix(61:62),lengthmatrix(61:62),'-*') 
hold  on 

polar(polarmatrix(63:64),lengthmatrix(63:64),':*r') 
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title(['16A{th}.  v_r=  ',num2str(doppler_fq(16)*c/fq_rr(l)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,17) 

polar(polarmatrix(65:66),lengthmatrix(65:66),'-*') 
hold  on 

polar(polarmatrix(67:68),lengthmatrix(67:68),':*r') 

title(['17A{th}.  v_r=  ',num2str(doppler_fq(17)*c/fq_rr(2)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,18) 

polar(polarmatrix(69:70),lengthmatrix(69:70),'-*') 
hold  on 

polar(polarmatrix(71:72),lengthmatrix(71:72),':*r') 

title(['18A{th}.  v_r=  ’,num2str(doppler_fq(18)*c/fq_rr(3)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,19) 

polar(polarmatrix(73:74),lengthmatrix(73:74),'-*') 
hold  on 

polar(polarmatrix(75:76),lengthmatrix(75:76),':*r') 

title(['19A{th}.  v_r=  ',num2str(doppler_fq(19)*c/fq_rr(l)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,20) 

polar(polarmatrix(77:78),lengthmatrix(77:78),'-*') 
hold  on 

polar(polarmatrix(79:80),lengthmatrix(79:80),':*r') 

title(['20A{th}.  v_r=  ',num2str(doppler_fq(20)*c/fq_rr(2)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,21) 

polar(polarmatrix(81:82),lengthmatrix(81:82),'-*') 
hold  on 

polar(polarmatrix(83:84),lengthmatrix(83:84),':*r') 

title(['21A  {th}.  v_r=  ',num2str(doppler_fq(21)*c/fq_rr(3)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,22) 

polar(polarmatrix(85:86),lengthmatrix(85:86),'-*') 
hold  on 

polar(polarmatrix(87:88),lengthmatrix(87:88),':*r') 

title(['22A{th}.  v_r=  ',num2str(doppler_fq(22)*c/fq_rr(l)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,23) 

polar(polarmatrix(89:90),lengthmatrix(89:90),'-*') 
hold  on 

polar(polarmatrix(91:92),lengthmatrix(91:92),':*r') 

title(['23A{th}.  v_r=  ',num2str(doppler_fq(23)*c/fq_rr(2)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,24) 

polar(polarmatrix(93:94),lengthmatrix(93:94),'-*') 
hold  on 

polar(polarmatrix(95:96),lengthmatrix(95:96),':*r') 
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title(['24A{th}.  v_r=  ',num2str(doppler_fq(24)*c/fq_rr(3)/2),'  m/s']) 

hline  =  findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,25) 

polar(polarmatrix(97:98),lengthmatrix(97:98),'-*') 
hold  on 

po  lar(po  larmatrix) 9 9 : 1 0 0 ) , 1 engthmatrix(9 9 : 1 0 0 ) , ' :  *  r') 

title(['25A{th}.  v_r=  ',num2str(doppler_fq(25)*c/fq_rr(l)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,26) 

polar(polarmatrix(  101:1 02  ),lengthmatrix(  101:1 02),'-*') 
hold  on 

polar(polarmatrix(  103 : 1 04),lengthmatrix(  103 : 1 04),':  *r') 

title(['26A{th}.  v_r=  ',num2str(doppler_fq(26)*c/fq_rr(2)/2),'  m/s']) 

hline  =  findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,27) 

polar(polarmatrix(  105:1 06),lengthmatrix(  105:1 06),'-*') 
hold  on 

polar(polarmatrix(  107:1 08), lengthmatrix(  107:1 08), ':*r') 

title(['27A{th}.  v_r=  ',num2str(doppler_fq(27)*c/fq_rr(3)/2),'  m/s']) 

hline  =  findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,28) 

polar(polarmatrix(  109:1 1 0),lengthmatrix(  1 09 : 1 1 0),'-*') 
hold  on 

polar(polarmatrix(  1 1 1:1 12),lengthmatrix(l  1 1 : 1 12),':*r') 

title(['28A{th}.  v_r=  ',num2str(doppler_fq(28)*c/fq_rr(l)/2),'  m/s']) 

hline  =  findobj(gca, 'Type', 'line'); 

set(hline,'LineWidth',3) 

subplot(5,6,29) 

polar(polarmatrix(  113:1 1 4),lengthmatrix(  113:114),'-*') 
hold  on 

polar(polarmatrix(  115:1 1 6),lengthmatrix(  115:11 6),':  *r') 

title(['29A{th}.  v_r=  ',num2str(doppler_fq(29)*c/fq_rr(2)/2),'  m/s']) 

hline  =  findobj(gca,'Type','line'); 

set(hline,'LineWidth',3) 

subplot(5,6,30) 

polar(polarmatrix(  117:11 8),lengthmatrix(  117:118),'-*') 
hold  on 

polar(polarmatrix(  119:1 20),lengthmatrix(  119:1 20),':  *r') 
title(['30A{th}.  v_r=  ',num2str(doppler_fq(30)*c/fq_rr(3)/2),'  m/s']) 
hline  =  findobj(gca,'Type','line'); 
set(hline,'LineWidth',3) 


mean(dopplerfq) 

std(dopplerfq) 

4)  WR_pulse_com 

%  Weather  Radar  Signal  Simulator,  developed  version 
%  Pulse  compression 
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%  Ulf  Schroder 
%  August,  2005 

% - 

%  Main  program  for  The  NPS/Marielle  Gosset  Weather  Radar  Signal  Simulator. 

% - 

%  DESCRIPTION 
%  Part  I 

%  From  Radar  Data: 

%  -  Generate  Beam  Resolution  cell  approximation 

%  -  Generate  Cube  adding  a  margine  to  the  Beam  Resolution  cell  to  allow  for 
%  the  rain  drops  to  fall  trough  the  through  beam  over  time.  The  stretched 
%  pulse  will  make  the  resolution  cell  larger  and  divided  into  range  bins 
%  according  to  specified  input  data.  To  get  a  resonable  average  the  middle 
%  range  bin  will  serve  as  the  "average"  range  bin 
% 

%  Part  II 

%  From  Rain  Parameters: 

%  Generating  the  fundamental  T  Matrix  concisting  of 
%  -  Diameter  Vector  and  Number  of  drops  in  each  diameter  interval 
%  using  the  Marshall-Palmer  exponential  approximation. 

%  -  the  backscatter  RCSs  for  each  Diameter  size 
%  -  the  total  number  of  drops  for  each  size  per  mA3 

%  Part  III 

%  Randomly  place  all  drops  in  cube 
%  Part  IV 

%  For  every  pulse  the  coherent  Electric  field  return  is  calculated.  Between 
%  pulses  all  drops  are  moved. 

%  Part  V 

%  Calculate  and  plot  results 

% - 

clear  all 
close  all 
clc 

%  Input  Data 

c=3e8;  %  Speed  of  light  [m/s] 

%  Radar  Data 

P_t=  1 0e3 ;  %  T ransmitted  power 

G_dB=20;  %  Maximum  gain  in  dB 

G=10A(G_dB/10);  %  Maximum  gain 

angle_3db=0.20;  %  3  dB  Beam  angle  [degrees],  (half  half  power  beamwidth) 

pulsewidth=0.20;  %  Pulsewidth  [microsecond] 

elevation=3;  %  Elevation  angle  [degrees] 

margine=5;  %  margin  for  resolution  cell[m] 

azimuth=0;  %  Azimuth  angle  [rad] 

PRT_l=l/200;  %  Pulse  Repetition  Time  [seconds] 

distance=0.5;  %  Range  from  the  radar  to  the  resolution  cell  [Km] 
wavelength=0.1;  %  Wavelength  [m] 

frequency=c/wavelength;  %  Frequency  [Hz] 
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% - 

%Introduced  in  Pulse  Compression  version 

chirp_width=2e6;  %  Frequency  step  in  [Hz] 

n_steps=l  0;  %  Number  of  Frequency  steps 

fq_step=chirp_width/n_steps;  %  Frequency  step  in  [Hz] 

d_range=c*pulsewidth*n_steps*le-6/2/n_steps;  %  range  bin  due  to  pulse  compression 

num_rans=10;  %  Number  of  runs  (pulse  compressed) 

range_bin=(distance*1000):d_range:(distance*1000+(n_steps-l)*d_range); 

% - 

%  Rain  parameters 

diam_low=0.05;  %  Lowest  diameter  of  Precipitation  [mm] 
diam_lim=l;  %  Divider  between  low  and  mid  [mm] 

diam_mid=5;  %  Divider  between  mid  and  high  [mm] 

diam_high=15;  %  Highest  diameter  of  Precipitation  [mm] 
resolutionl=0.08;  %  Lower  diameter  resolution  [mm] 
resolution2=0.3;  %  Higher  diameter  resolution  [mm] 
resolution3=0.8;  %  Highest  diameter  resolution  [mm] 
sort_approx=l;  %  RCS  approximation  model;  [l]=Rayleigh 
kw=0.93;  %  abs(K)A2  approx  for  water 

rain_rate=50;  %  Rain  Rate  [mm/h] 

rain_anglel=0;  %  degrees  Rain  Angle  (Phi) 

rain_angle2=0;  %  degrees  Rain  Angle  (Theta) 
wind_vector=[-10;0;0];  %  Wind  speed  vector  (x,y,z)  [m/s] 
spread=l;  %  Wind  speed  spread  [m/s],  (all  directions) 

spread_2=4;  %  Drop  velocity  spread  [m/s] 

Norm=200;  %  Normalization  coefficient 

%  Directional  cosines 

u=sin(pi/2-elevation*pi/l  80)*cos(azimuth*pi/l  80); 
v=sin(pi/2-elevation*pi/180)*sin(azimuth*pi/180); 
w=cos(pi/2-elevation*pi/180); 

% - 

%  Part  I,  Generate  the  resolution  volume 

% - 

%  Approximation  of  the  resolution  cell  of  the  pulse  circular  box, 

%  for  any  azimuth  angle 

% - 

%  INPUT 
%  -  distance 
%  -  angle_3dB 
%  -  pulsewidth 
%  -  elevation 
%  -  margine 
%  -  azimuth 

% - 

%  OUTPUT 

%  -  volume  resolution:  computed  volume  of  the  resolution  cell 
%  -  volume_box:  volume  of  the  box  where  the  drops  are  put 
%  -  box:  the  8  rectangular  coordinates  of  the  box 
%  -  coord_vol_res:  the  coordinates  of  the  resolution  cell  where  the  drops 
%  are  taken. 

% - 
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% - 

%  Creating  the  circular  beam  resolution  cell 
range_res=(c*pulsewidth*n_steps*le-6)/2; 


%  small  angle  approximation  allows  for  following 
angle_rad=(angle_3  db  *pi)/ 1 8  0 ; 

%  3  dB  beamwidth  [rad] 
radius  l=(tan(angle_rad)*(distance*  1000)); 

%  Radius  of  beam  at  R=distance 
radius2=(tan(angle_rad)*((distance*  1000)+range_res));  ... 

%  Radius  of  beam  at  R=distance  +  range  resolution 

surface  1  =pi  *radius  1 A2 ; 
surface2=pi*radius2A2; 

%  volume  of  cone  =  base  surface  times  height  for  all  three 
volume_conel=(surfacel*(distance*1000))/3; 
volume_cone2=(surface2*((distance*1000)+range_res))/3; 
volume_resolution=(volume_cone2-volume_conel);  %  Beam  Cone  Volume 

%  To  assure  right  volume  for  averaging 
DR=1:10; 

range_res_ref=(c*pulsewidth*DR*  le-6)/2; 

radius2_ref=(tan(angle_rad)*((distance*1000)+range_res_ref));  ... 
surface2_ref=pi*radius2_ref.A2; 

volume_cone2_ref=(surface2_ref.*((distance*1000)+range_res_ref))/3; 

volume_resolution_ref=(volume_cone2_ref-[volume_conel  volume_cone2_ref(l:end-l)]); 

%  Beam  Cone  Volume 


%  computation  of  the  resolution  cell  coordinates,  spherical  coordinates 
%  elevation  in  radians 

elevation_rad=(elevation*pi)/180;  %  Elevation  [rad] 

azimuth_rad=(azimuth*pi)/180;  %  Azimuth  [rad] 

%  Matrix  used  when  comparing  wheather  a  drop  is  inside  or  outside 
%  beamwidth 

coord_vol_res=[distance*  1000,  (distance*  1000)+range_res;elevation_rad-angle_rad,... 

elevation_rad+angle_rad;azimuth_rad-angle_rad,azimuth_rad+angle_rad]; 

% -  - - 


% - 

%  Creating  the  box 

%  computation  of  the  box  and  its  coordinates 

coin  1  =[azimuth_rad+angle_rad,elevation_rad+angle_rad, distance*  1000]; 
coin2=[azimuth_rad-angle_rad,elevation_rad+angle_rad,distance*  1 000] ; 
coin3=[azimuth_rad+angle_rad,elevation_rad-angle_rad,distance*1000]; 
coin4=[azimuth_rad-angle_rad,elevation_rad-angle_rad,distance*1000]; 
coin5=[azimuth_rad+angle_rad,elevation_rad+angle_rad,distance*1000+range_res]; 
coin6=[azimuth_rad-angle_rad,elevation_rad+angle_rad,  distance*  1000+range_res]; 
coin7=[azimuth_rad+angle_rad,elevation_rad-angle_rad, distance*  1000+rangeres]; 
coin8=[azimuth_rad-angle_rad,elevation_rad-angle_rad,distance*1000+range_res]; 

[x  1  ,y  1  ,zl  ]=sph2cart(coin  1  ( 1  ),coin  1  (2), coin  1  (3 )); 
[x2,y2,z2]=sph2cart(coin2(l),coin2(2),coin2(3)); 
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[x3,y3,z3]=sph2cart(coin3(l),coin3(2),coin3(3)); 

[x4,y4,z4]=sph2cart(coin4(l),coin4(2),coin4(3)); 

[x5,y5,z5]=sph2cart(coin5(l),coin5(2),coin5(3)); 

[x6,y6,z6]=sph2cart(coin6(l),coin6(2),coin6(3)); 

[x7,y7,z7]=sph2cart(coin7(l),coin7(2),coin7(3)); 

[x8,y8,z8]=sph2cart(coin8(l),coin8(2),coin8(3)); 

xmax=max([xl,x2,x3,x4,x5,x6,x7,x8]); 
xmin=min([xl,x2,x3,x4,x5,x6,x7,x8]); 
ymax=max([y  1  ,y2,y3,y4,y5,y6,y7  ,y8]); 
ymin=min([y  1  ,y2,y3,y4,y5,y6,y7  ,y8]); 
zmax=max([zl,z2,z3,z4,z5,z6,z7,z8]); 
zmin=min([zl,z2,z3,z4,z5,z6,z7,z8]); 

%  the  box  with  its  margin 

coin_  1  =[xmin-margine,ymin-margine,zmin-margine] ; 

coin_2=[xmin-margine,ymax+margine,zmin-margine]; 

coin_3=[xmax+margine,ymax+margine,zmin-margine]; 

coin_4=[xmax+margine,ymin-margine,zmin-margine]; 

coin_5=[xmin-margine,ymin-margine,zmax+margine]; 

coin_6=[xmin-margine,ymax+margine,zmax+margine]; 

coin_7=[xmax+margine,ymax+margine,zmax+margine]; 

coin_8=[xmax+margine,ymin-margine,zmax+margine]; 

box=[coin_r,coin_2',coin_3',coin_4',coin_5',coin_6',coin_7',coin_8']; 

%  calculate  the  volume  of  the  box 

volume_box=abs((xmax+margine-(xmin-margine))*(ymax+margine-(ymin-margine))*... 

(zmax+margine-(zmin-margine))); 


cube=[coin_r,coin_2',coin_3',coin_4',coin_r,coin_5',coin_6',coin_2',... 

coin_6',coin_7',coin_3',coin_7',coin_8',coin_4',coin_8',coin_5']; 

% - 

% - 

%  Part  II,  Generate  the  fundamental  T  Matrix 

% - 

%  This  part  produces  a  matrix  of  all  drops  that  contribute  to  the  total 
%  reflectivity  of  the  resolution  cell 

% - 

%  First  all  sizes  are  created,  with  two  vectors  of  different  size  resolution. 

%  The  size  vectors  represents  all  diameters.  Next  the  backscatter  RCS  for  each 
%  size  is  computed. 

%  Using  the  Marshall-Palmer  exponential  approximation,  the  number  of  drops  of 
%  each  diameter  is  computed. 

% - 

%  INPUT 
%  -  diam_low 
%  -  diam_lim 
%  -  diam_high 
%  -  resolution  1 
%  -  resolution2 
%  -  sort_approx 
%  -  wavelength 


142 


%  -  Kw 
%  -  rain_rate 

% - 

%  OUTPUT 
%  -  number_drops 
%  Matrix  T : 

%  -  first  row  all  sizes 
%  -  second  row  all  the  backscatter  RCSs 
%  -  third  row  the  total  number  of  drops  for  each  size  per  mA3 
%  the  columns  start  with  the  smallest  size 

% - 

%  Create  dropsize  vector 

%  Precipitaion  Diameter  vector  (lowest<D<divider) 
vector_size_l=diam_low:resolutionl:diam_lim; 

%  Precipitaion  Diameter  vector  (divider<D<mid) 
vector_size_2=diam_lim:resolution2:diam_mid; 

%  Precipitaion  Diameter  vector  (mid<D<highest) 
vector_size_3=diam_mid:resolution3:diam_high; 

%  Precipitaion  Diameter  Vector 

vector_size=[vector_size_l,vector_size_2,vector_size_3]; 


% - 

%  Create  RCS  matrix  for  all  Rain  drops  of  Precipitation  Vector  and  for  all 
%  frequencies. 

fq_rr= frequency  :fq_step:frequency+(n_steps-l)*fq_step; 
rcs(l:n_steps,:)=((kwA2)*(piA5)/(cA4)*fq_rr'.A4)*((vector_size*le-3).A6); 
% - 


% - 

%  The  Marshall  Palmer  approximation  is  used  to  create  Precipitaion  Diameter 
%  Distribution  Vector  (Drop  Size  Distribution  N(D)  Vector) 

%  The  number  of  drops  are  scaled  down  by  a  factor  (Norm)  to  improve  the 
%  speed  of  the  simulator 

% - 

lambda=4.1*(rain_rateA(-0.21));  %[mm-l] 

number=8000*exp(-lambda*vector_size)/Norm; 

vector_resl=ones(l,length(vector_size_l))*resolutionl; 
vector_res2=ones(l,length(vector_size_2))*resolution2; 
vcctor_rcs3~oncs(  1 .  length)  vcctor_sizc_3))*reso  I  ution3; 

%  delta(Diameter)  dD 

vector_res=[vector_resl,vector_res2,vector_res3]; 

%  the  approximate  "true"  number  of  drops  N(D)dD  vector 
number=number.*vector_res; 

% - 


%■ 


%  create  matrix  T  without  res 
T=[vector_size;number]; 


%■ 
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% - 

%  Change  matrix  T  to  matrix  'matrixTRound' 

% - 

%  This  part  computes  the  total  number  of  drops  by  utilizing  the  resolution 
%  cell  and  the  matrix_T_Round 

%  The  only  difference  between  the  matrix_T_Round  and  T,  is  the  fact  that 
%  matrix_T_Round  has  the  total  number  of  drops  on  the  third  row  and  they 
%  are  rounded  of  to  nearest  integer  value. 

% - 


matrix_T_Round=T ; 

matrix_T_Round(2,:)=round(matrix_T_Round(2,:)*volume_box); 

%  Keep  track  of  errors  due  to  round  off 
diff=T(2,:)*volume_box-matrix_T_Round(2,:); 

error_rcs=diff*rcs'*Norm/volume_box*sum(volume_resolution_ref)/n_steps; 

error_Z=error_rcs/sum(volume_resolution_ref)/n_steps*cA4./fq_rr.A4/piA5/... 

abs(kw)A2/le-18; 


% - 

%  Part  III,  Initial  positioning  of  the  drops  within  the  volume 

% - 

%  This  part  generates  the  initial  position  of  the  drops  in  the  volume  of 
%  the  box. 

%  The  positions  are  recorded  in  the  file  "sizexx.maf  taking  into  account 
%  the  azimuth  angle 

% - 

%  INPUT 
%  -  box 

%  -  matrix_T_Round 
%  -  elevation 
%  -  azimuth 

% - 

%  OUTPUT 

%  All  files  including  the  sizes.  Each  contains  the  positions  of  all  drops 
%  for  a  given  diameter.  The  data  recoreded  as  follows: 

%  -  radial  distance  r 
%  -  elevation  phi 
%  -  azimuth  theta 

%  -  angle  with  respect  to  beam  center.  This  parameter  will  be 
%  used  to  identify  drops  which  are  in  the  beam. 

% - 


%  Randomly  put  the  calculated  number  of  drops  of  each  diameter  in  the  box 
number_size=length(matrix_T_Round(l,:)); 

for  i=l:number_size 

%  select  number  of  drops  of  certain  Diameter 
drop  s=matrix_T_Round(2 ,  i ) ; 

%  create  as  many  positions  on  the  x  axis  as  #  drops  per  size 
elements_x=rand(  1  ,drops); 

%  put  them  in  the  range  of  interest 
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%  Random  placement  x  for  every  drop 

position_x=box(  1 , 1  )+(box(  1 ,3)-box(  1 , 1  ))*elements_x; 

%  create  as  many  positionson  the  y  axis  as  #  drops  per  size 
elements_y=rand(  1  ,drops); 

%  put  them  in  the  range  of  interest 
%  Random  placement  y  for  every  drop 

position_y=box(2, 1  )+(box(2,2)-box(2, 1  ))*elements_y; 

%  create  as  many  positions  on  the  z  axis  as  #  drops  per  size 
elements_z=rand(  1 ,  drops); 

%  put  them  in  the  range  of  interest 
%  Random  placement  z  for  every  drop 

position_z=box(3,l)+(box(3,5)-box(3,l))*elements_z; 

%  matrix  of  results 

%  Creating  a  drop  Matrix  using  spherical  coordinates 

[theta,  phi,r]=cart2sph(position_x,position_y,position_z); 

%  Calculating  the  range,  phi,  theta  and  angle  to  every  drop  with 
%  reference  to  phi_0  and  theta_0.  Used  to  estimate  weighted  power  return. 
vector_azimuth=ones(l,drops)*((azimuth*pi)/180); 
vector_elevation=ones(  1  ,drops)*((elevation*pi)/ 180); 
[xf,yf,zfj=sph2cart(vector_azimuth,vector_elevation,r); 
range=sqrt((position_x-xf).A2.+(position_y-yf).A2.+(position_z-zf).A2); 

%  approximation 
angle=atan(range./r); 

%  Saving  the  drop  positions  in  size%d.mat 
%  matrixpos 

matrixpos=[r;phi;theta;angle]; 
save(sprintf('size%d. mat'd), 'matrixpos'); 
end 

% - 


% - 

%  Part  IV,  Calculate  Coherent  Electric  field  return. 
% - 


%  For  every  pulse  the  coherent  Electric  field  return  is  calculatet.  Between 
%  pulses  all  drops  are  moved. 


% 


%  Run  the  simulation  for  specified  number  of  pulses 
% - 


%  Zero  all  vectors  and  matrices  that  are  going  to  be  used 
phase_vector=[]; 

sum_E=zeros(num_runs,n_steps); 

weighted_E=zeros(num_runs,n_steps); 

mnning_Z=[]; 

%  Initiate  reference  variables 
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%  Used  to  separate  inns 


p=0; 

for  runs=l  :num_runs 
p=p+l; 

hh=0;  %  Used  to  separate  frequencies 

for  fq_new=frequency:fq_step:frequency+(n_steps-l)*fq_step 

hh-hh+1 ; 

k=2*pi*fq_new/c; 

disp(sprintf('calculate  the  power-return  of  pulse  %d\n',p)); 


% - 

%  computes  the  (weighted)  I  and  Q  return 

% - 

%  computes  the  total  number  of  drops 

% - 

%  INPUT 
%  -  Matrix  T 
%  -  wavelength 
%  -  angle_3dB 
%  -  distance 
%  -  coord_vol_res 
%  -  number_pulses 

% - 

%  OUTPUT 

%  -  totalE:  the  received  complex  Electric  field) 

%  -  power:  the  power 
%  -  voltage:  the  voltage 

%  -  drops_beam:  the  number  of  drops  in  the  beam 

% - 

%  Only  the  drops  which  are  in  the  resolution  cell  are  selected  for  the 
%  computation 

% - 


number_size=length(matrix_T_Round(  1 ,:)); 

total_E=0; 

total_phase=0; 

new_E=0; 

for  i=l:number_size 

%  opening  the  file  of  registered  data  matrixpos. 
load(sprintf('size%d.mat',i)); 

%  removing  drops  that  are  outside  the  resolution  cell 

number_drops_per_size=length(matrixpos(l,:)); 

dropmatrix=[]; 

%  Filling  the  drop  matrix  with  all  drops  inside  the  beam 

%  resolution  volume 

for  d=l  :number_drops_per_size 

if  ((matrixpos(4,d)<=(angle_3db*pi/180))&(matrixpos(l,d)>=... 

range_bin(hh))&(matrixpos(l,d)<=i'ange_bin(hh)+d_range)); 
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dropmatrix=[dropmatrix,matrixpos( :  ,d)] ; 
end 
end 


% - 

%  Convert  to  cartesian  coordinates 
if  isempty(dropmatrix)==0 

[position_x,position_y,position_z]=sph2cart(dropmatrix(3 . 
dropmatrix(2,  :),dropmatrix(  1, :)); 

%  Create  comparison  vectors 

vector_azimuth=ones(l,length(dropmatrix(l,:)))*((azimuth*pi)/180); 

vector_elevation=ones(l,length(dropmatrix(l,:)))*((elevation*pi)/180); 


%  Create  a  phase  vector  for  the  drops  of  current  Diameter 

r_phase=position_x.*u+position_y.*v+position_z.*w; 

phase=exp((2*j*k).*r_phase); 

%  Sum  the  contributions  to  the  E-field 

total_phase=total_phase+sum(phase); 

total_E=total_E+sum(sqrt(rcs(hh,i)).*phase); 


% - 

%  ceate  weighted  power  return,  concidering  range  (r),  theta,  phi 
%  with  reference  to  boresight  and  radar  parameters 

E_weighted=sqrt(P_t)*G*c/fq_rr(hh)/(sqrt(4*pi)A3)./dropmatrix(l,:).A2.*... 

(exp(-4*log(2).*((dropmatrix(3,:)-vector_azimuth).A2/... 

(2*angle_3db*pi/180)A2+((dropmatrix(2,:)- 

vector_elevation).A2/(2*angle_3db*pi/180)A2)))); 

new_E=new_E+sum(E_weighted.*(sqrt(rcs(hh,i)).*phase)); 

end 

end 


% - 

%  add  to  E-field  and  phase  vector 
sum_E(p,hh)=total_E; 
weighted_E(p,hh)=new_E; 
phase_vector(2*p-l)=ANGLE(total_phase); 
% - 


end 

% - 

%  This  part  simulates  the  movement  of  the  drops  and  gives  new  positions. 
%  First  using  the  Doppler  PRT  (PRT1) 

% - 

%  INPUT 

%  -  matrix  T  Round 
%  -  PRT 
%  -  azimuth 
%  -  elevation 
%  -  rain_anglel 
%  -  rain_angle2 
%  -  diam_lim 
%  -  resolution! 
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%  -  resolution2 

%  -  windvector 

% - 

%  OUTPUTS 

%  Output  files  containing  sizes  with  new  positions 
%  structured  the  same  way  as  in  previous 

%  there  is  also  an  output  file  for  speed  (one  for  each  size)  for  the  pulse 
%  pair. 

% - 


number_size=length(matrix_T_Round(l,:)); 

%  adding  a  Gaussian  random  term  to  make  the  wind  speed  spread 
if  wind_vector(l)~=0 

new_wind_vector(  1  )=wind_vector(  1  )*normrnd(  1  ,abs(spread/max(wind_vector(  1 )))); 
else 

new_wind_vector(  1  )=0; 
end 

if  wind_vector(2)~=0 

new_wind_vector(2)=wind_vector(2)*normrnd(l,abs(spread/max(wind_vector(2)))); 

else 

new_wind_vector(2)=0; 

end 

if  wind_vector(3)~=0 

new_wind_vector(3)=wind_vector(3)*normrnd(l,abs(spread/max(wind_vector(3)))); 

else 

new_wind_vector(  3  )=0; 
end 

for  i=T:number_size 

%  opening  of  data  file  registered  in  matrixpos. 
load(sprintf('size%d.mat',i)); 

%  total  number  of  drops 
drops=matrix_T_Round(2,i); 

%  transforming  to  rectangular  coordinates 
[x,y,z]=sph2cart(matrixpos(3,:),matrixpos(2,:),matrixpos(l,:)); 

diameter=matrix_T_Round(  1  ,i); 


% - 

%  Drop  fall  speed  depends  on  drop  diameter  in  the  z  axis, 

%  Atlas-Ulbrich  approximation  is  used 

%  Approximation  is  valid  in  the  diameter  range  5*e-4,  5*e-3 
% - 


wtmax=386.6*((matrix_T_Round(l,i)*0.001)A0.67); 

velocityz=(wtmax*(ones(l,drops)))*cos((rain_anglel*pi)/180); 

velocityx=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 

cos((rain_angle2*pi)/180); 

velocityy=(wtmax*(ones(l,drops)))*sin((rain_anglel*pi)/180)*... 
sin((rain_angle2*pi)/l  80); 
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% - 

% - 

% 

%  %  adding  another  Gaussian  random  term  to  make  the  wind  speed 

%  %  spread  between  different  drops  sizes 

%  if  new_wind_vector(  1  )~=0 

%  newest_wind_vector(  1  )=new_wind_vector(  1  )*normrnd(  1  ,abs(spread. 

%  /max(wind_vector(l)))); 

%  else 

%  newest_wind_vector(  1  )=0; 

%  end 

%  if  new_wind_vector(2)~=0 

%  newest_wind_vector(2)=new_wind_vector(2)*normrnd(l,abs(spread. 

%  /max(wind_vector(2)))); 

%  else 

%  newest_wind_vector(2)=0; 

%  end 

%  if  new_wind_vector(3)~=0 

%  newest_wind_vector(3)=new_wind_vector(3)*normrnd(l,abs(spread. 

%  /max(wind_vector(3)))); 

%  else 

%  newest_wind_vector(3)=0; 

%  end 

% 

%  %  If  spread  within  pulse 

%  velocityxt=newest_wind_vector(  1  )-velocityx; 

%  velocityyt=newest_wind_vector(2)-velocityy; 

%  velocityzt=newest_wind_vector(3)-velocityz; 

% 

% 

% - 

% - 

%  If  spread  for  all  drops  regardless  of  size 
newest_wind_vector=[] ; 


%  adding  Gaussian  random  term  to  make  the  wind  speed 
%  spread  between  drops 
if  new_wind_vector(l)~=0 

ne  west_wind_vector(  1 , :  )=new_wind_vector(  1 )  * . . . 

normrnd(l,abs(spread_2/max(new_wind_vector(l))),l,  drops); 


else 

newest_wind_vector(l,:)=0; 

end 

if  new_wind_vector(2)~=0 

newest_wind_vector(2,:)=new_wind_vector(2)*... 

normrnd(  1  ,abs(spread_2/max(new_new_wind_vector(2 ))),  1  ,drops); 


else 

newest_wind_vector(2,:)=0; 

end 

if  new_wind_vector(3)~=0 
ne  west_wind_vector(3 , :  )=new_wind_vector(  3 )  * . . . 

normrnd(l,abs(spread_2/max(new_wind_vector(3 ))),!,  drops); 


else 

ne  west_wind_vector(  3 , :  )=0; 
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end 


%  If  spread  within  pulse 
velocityxt=newest_wind_vector(l,:)-velocityx; 
velocityyt=newest_wind_vector(2,:)-velocityy; 
velocityzt=newest_wind_vector(3 , :  )-velocityz; 


%■ 

%• 


%  speed  registration 

m_vclocity— |  velocityxt;vclocityyt;velocityzt  |; 
save(sprintf('velocity%d.mat',i),'m_velocity'); 

deplacementx=velocityxt*PRT_l ; 
deplacementy=velocityyt*  PRT_  1 ; 
deplacementz=velocityzt*PRT_l ; 

position_x=x+deplacementx; 

position_y=y+deplacementy; 

position_z=z+deplacementz; 

%  transformation  into  spherical  coordinates 

[theta,  phi,r]=cart2sph(position_x,position_y,position_z); 

%  computation  of  new  angles 
vector_azimuth=ones(l,drops)*(azimuth*pi)/180; 
vector_elevation=ones(l,drops)*(elevation*pi)/180; 
[xf,yf,zf]=sph2cart(vector_azimuth,vector_elevation,r); 

range=sqrt((position_x-xf).A2.+(position_y-yf).A2.+(position_z-zf).A2); 
%  approximation 
angle=atan(range ./ r) ; 

%  matrixpos 

matrixpos=[r;phi;theta;angle]; 

%  save  to  file,  matrixpos 

s=sprintf(  'size%d.  mat' ,  i) ; 
save(s,  'matrixpos'); 


end 

end 


%• 


% - 

%  calculate  RCS,  Z  and  Power  return  for  every  pulse  including  error.  The 
%  reflectivity  is  caculated  for  each  volume  bin  and  then  summed  and 
%  averaged. 

% - 

for  o=l  :num_runs 

running_Z(o,:)=(Norm*abs(sum_E(o,:)).A2)*cA4./fq_rr(o)A4/piA5/... 

abs(kw)A2/volume_resolution_ref(o)/le-18+error_Z(o); 

end 
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running_Z=sum(running_Z,2)7n_steps; 


%  Transform  into  a  vector 


% - 

%  Adding  the  error  to  the  weighted_E  estimations 

% - 

error_power=error_rcs*P_t*GA2*cA2./fq_rr.A2/piA2*(2*angle_3db*pi/180)A2*... 

(c*pulsewidth*le-6)/1024/log(2)/piA2/(distance*le3)A2; 

for  o=l  :num_mns 

power_return(o,:)=Norm*abs(weighted_E(o,:)).A2+error_power(o); 

end 

power_retiirn=sum(power_return,2)7n_steps;  %  Transform  into  a  vector 

%  to  use  to  estimate  zeroth  moment 
for  o=l:n_steps 

power_return_ref_matrix(:,o)=Norm*abs(weighted_E(:,o)).A2+error_power(o); 

end 

power_return_ref=sum(power_return_ref_matrix,l)/num_runs;  %  Transform  into  a  vector 


% - 

%  calculate  average  Power  and  Z 

% - 

Z_avg=sum(running_Z)/p; 

Z_avg_plot=ones(  1  ,p)*Z_avg; 

power_return_avg=sum(power_return)/p;  %  [W] 

power_return_avg_plot=ones(l,p)*power_return_avg; 

number_vector=l  :p; 


% - 

%  Calculate  references 
%  Integral  form  of  Z 

% - 

ii=0: 100/1000: 100; 

zz=ii.A6.*8000.*exp(-4.1.*rain_rate.A(-.21).*ii); 

real_Z=ones(l,p)*10*logl0(trapz(ii,zz)); 


% - 

%  Evaluate  average  Power  return  to  make  Z  estimation  including  error 

% - 

%  To  make  the  estimate  with  the  correct  frequency,  the  power  return  Matrix 
%  is  used.  For  the  plot  the  result  must  me  summed  and  divided  by  the  number 
%  of  samples  used  to  average  (which  means  over  both  frequency  and  pulse). 

Z_estimation=power_retum_ref_matrix*((fq_rr.A2).*((c./fq_rr).A4))7P_t/GA2/cA2*... 

piA2/(2*angle_3db*pi/180)A2/(c*pulsewidth*le- 

6)*1024*log(2)*(distance*le3+range_res/2)A2/... 

((kwA2)*(piA5))/le-18; 

Z_estimation_plot=ones(l,p)*10*logl0(sum(Z_estimation)/n_steps/num_runs); 


%■ 
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%  Part  V,  Plot  of  results 
% - 


figure(l); 

kl=l:p; 

number_vector=kl.*PRT_l  *  1 000; 

text_string=['Power  Return  for  all  ',num2str(p),'  pulse  trains.'... 

Frequency  stepped  from  ',num2str(fq_rr(l)/le9),'  GFlz  to  ',... 
num2str(fq_rr(n_steps)/le9),'  GFlz  with  \num2str(fq_step/le6)... 

MFlz  steps.']; 

%  also  converting  to  dBm 

plot(number_vector,  1 0*log  1 0(power_return)+30,'-xb',  number_vector,  10*logl 0.. . 

(power_retum_avg_plot)+30,'~  r'); 
title(text_string); 
xlabel('Sample  time  [ms]'); 
ylabel('Power  [dBm]'); 

legend('Power  return',  'Average  Power  return', 0) 
grid  on 

figure  (2) 

plot(number_vector,  1 0*log  1 0(abs(running_Z)),'-xb',number_vector,  10*... 
loglO(abs(Z_avg_plot)),'— r',number_vector,real_Z,'— k',... 
number_vector,Z_estimation_plot,'~ b'); 
title('Estimated  Z  per  pulse'); 
xlabel('Sample  time'); 
ylabel('Z  [dBZ]'); 

legend('Simulation  Z,  (instantaneous)',  'Estimated  Z,  (average  return)', 'Z=\int  DA6N(D)dD',... 

'Estimated  Z,  Zeroth  moment', 0) 
grid  on 
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